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I. INTRODUCTION 


1.1 Introduction 

Among many geometrical profiles which generate hydrodynamic pressure 
in fluid- film bearings, the step geometry is one of the most proficient 
methods to achieve this purpose. This discovery was first made by Lord 
Rayleigh in 1918 (Ref. 1 ), and subsequently there have been many 

contributions (Refs. 2 to 8 ) on the prediction of the performance of 

the Rayleigh-step thrust as well as journal bearings. 

Recently, there has been another significant application of the 
Rayleigh-step bearing in the field of dynamic sealing in advanced air- 
breathing propulsion systems (Ref. 9 ). In this application, a series 

of Rayleigh-step pads is employed on the high-pressure side of a face 
seal in order to maintain a small steady gap (in the order of 0.001") 
between a one-tooth laybrinth and the high- speech rotor surface (Fig. 1 ), 
The flat step is shrouded in order to minimize the side leakage. These 
pads function strictly as hydrodynamic thrust bearings operate in the 
high-ambient pressure , and provide the necessary stiffness to maintain a 
steady gap. The satisfactory operation of this type of seals depends 
critically on the dynamic performance of these thrust pads in the 
presence of an oscillatory force or a disturbance due to rotor run-out 
or rotor unevenness. 

So far, the investigations on shrouded, Rayleigh-step bearings have 
been restricted to the prediction of the static performance (Ref. 10 ) 
only. The dynamic characteristics of this type of thrust-pads have not 
been given much attention. It is the main objective of this report to 
study the influence of the nonlinear, gas- film, restoring and damping force 
upon the response of the pad to a given forcing function or disturbance. 


1 
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Specifically, the work reported herein fulfills the following 
objectives : 

1. To develop a gas- film analysis of purely hydrodynamic, 

Rayleigh-s tep pad to calculate the quasistatic stiffness and 
damping, which depend not only on the operating conditions 
but also on the vibration of the system. 

2. To extend the analysis to the nonlinear axial response of 
the stationary ring due to any external force excitation or 
any disturbances induced by the rotor misalignment or surface 
waviness. 

3. To develop a stability analysis of the infinitely 

wide, single-step pad to explore whether there exists any 
thresholds of stability for the system. 

1 . 2 Historical Survey 

In 1918, after Reynolds' theory of thick film lubrication became 
generally accepted. Lord Rayleigh (Ref. 1 ) first applied the theory and 
discovered that an optimum profile for the load capacity of a slider 
bearing is a flat step. The method of the calculus of variations was 
used to optimize the shape of slider bearing to an infinitely long, 
incompressible film. From then on, any hydrodynanic bearings composing 
of two sections of parallel surface film are called Rayleigh- step bearings. 

The optimized geometry of a Rayleigh-step slider and its corresponding 
optimized, non-dimensional, load capacity for an incompressible film can 
be found in most textbooks on lubrication (Ref. 11 ). For compressible 
gas film, the optimized, one-dimensional, Rayleigh- step bearing was analyzed 
by Wylie and Maday in 1970 (Ref. 2). The load capacity of an optimized step 



bearing was found to be slightly lower at low bearing number but much 
higher at higher bearing number than the load capacity of an incompressible, 
optimize^ slider bearing. 

The problem of a two-dimensional, Rayleigh- step, thrust bearing has 
not received much attention until 1954, C. F. Kett leborough (Ref. 3) 
solved the pressure profile and calculated the load capacity for the step 
thrust bearing by Relaxation methods. In 1955 (Ref.f7)he applied the 
analogy method contributed by A. Kingsbury (Ref. 12) to investigate the 
pressure profile for an oil-lubricated step bearing by an electrolytic tank. 

In 1959, using air as the lubricant, K. C, Kochi (Ref. 6) showed the 
characteristics of an in finitely- wide, Rayleigh- step, thrust pad by the 
use of semi-graphical method. He demonstrated that an analytical solution 
to express the pressure profile explicitly is extremely difficult, because 
the Reynolds equation for a compressible film is nonlinear. 

In 1961, J. S. Ausmann (Ref. 7) made certain approximations to 
linearize the Reynolds equation for a compressible film, and obtained a 
series solution of the pressure and load for a self-acting, stepped, sector, 
thrust bearing by the aid of Neumann polynomial. He obtained the numerical 
solutions for the optimized number of sectors and step depth for the maxi- 
mum load carrying capacity. 

Recently, Cheng, Chow and Wilcock (Ref. 8) presented some results 
for shrouded, Rayleigh- step pad used as a flexible face seal. The pres- 
sure generation and static stability of this type of surface profiles 
using as a flexible seal ring was discussed, and the effectiveness of 
hydrodynamic action was confined to the static stiffness characteristics 
of gas film. The influence of the nonlinear gas- film forces and the 
question of gas-film stability of the hydrodynamic, Rayleigh-step pad 
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as not been investigated. Therefore, to ensure a safe operation of 
ayleigh-step, thrust pads, it ia necessary to conduct a full-scale nonlinear 
tudy, which is the major objective of this work. 



I I • GAS FILM FORCE S 


2.1 Statement of the Problem 

The major concern in this section is the determination of the 
time- dependent, pressure distribution within the gas film between two 
annular surfaces containing a series of Rayleigh-step pockets as 
shown in Fig. 1. This problem is formulated within the framework of 
the conventional theory of lubrication for a compressible lubricant. 

The major assumptions commonly used for gas-lubrication theories are: 

1. The pressure across the film is constant 

2. The flow is laminar 

3. Inertia forces are neglected 

4. The film is isothermal 

5. The flow is Newtonian. 

Under these assumptions, the governing equation for a transient, 
continuous, pressure field becomes the well-known, transient, Reynolds 
equation for a compressible fluid. Various analytical and numerical 
methods for the solution to this equation have been outlined by 
Castelli and Peviric (Ref. 13 ). For the present problem, the abrupt 

geometry introduced by the Rayleigh pocket makes it rather difficult to 
solve the pressure equation by any analytical methods. For this reason, 
a numerical method is employed in solving the discretized pressure 
field based on the solution of a system, non-linear algebraic equations 
by the Newton-Raphson procedure. The numerical integration of the 
discretized pressure field gives rise to the time -dependent gas-film 
forces which are required for investigating the nonlinear response 
of this type gas-film to a prescribed forcing function or disturbance 
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function. 




2.2 Pressure Equation 

The discretized pressure equation is derived by considering a flow 
balance within an element of the gas film as shown in Fig. 2. The mass 
flow rate into the left boundary, AD, and the bottom boundary, AB, are 
designated as and q 2 * The mass flow rates out of the boundaries, 

DC and BC are, likewise, designated as q^ and q^. In addition, there are 
mass stored in the volume which is designated as q^. 

Based on the assumption in 2.1, the velocity Is parabolic across 
the film. Integrating the velocity across the film, one may express 
the mass flow rate in terms of the boundary velocities and the pressure 
gradients. These can be written as 


[p JlEh . iE 
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r U h .3 
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The flow balance requires 
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Introducing the following nondimens tonal variables: 
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and using the equation of state p * , Equation (2.6) becomes 
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Further simplification leads to. 
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Introducing the following finite difference approximations 
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one obtains a system of first order differential equations in time for 
the discretized pressures. 

To obtain the exact, time-dependent, gas-film forces, it is necessary 
to solve this equation simultaneously with the equations of motion of 
the supported mass by an explicit or implicit, numerical procedure for the 
initial-valued problems. This procedure necessitates the calculation of 
the pressure field at each time interval during the transient, and it is 
extremely uneconomical and cumbersome. 

In the case of a high, ambient pressure, the effect produced by the 

pip iu 

term containing ~ becomes insignificant comparing to that by — , one may 
ft? 

neglect the — effect, and thus decoupled the gas- film force calculation 
from the dynamics of the supporting mass. Ignoring the contribution by 

~ , the pressure distribution can be solved independently as a function 

ftH * ftp 

of H and ~ or H. The assumption of a negligible " effect is made in the 

subsequent analysis. Expanding the terms in Equation (2.7) by using 

relations (2.8) one obtains 


l 3 Q i,j-l + a l A.J-l + a 5 Q i-l,j ■ (a 3 + a 4 + a 5 + a 6 )Q i,j + a 6 Q i+l,j 
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2.3 Method of Solution 

The set of nonlinear, algebraic equations to be solved for the pressure 
field, Q , is 


'i,j * a 3 Q i,j-l + a l A.j-1 + 


■sVl,] " (a 3 + a 4 + a 5 + a 6 )Q i,j + LH~~ “ 

1 * J “ 


+ a 6 Q i+l , j + a 4 Q l,j+l ' a 2^1,j+l = 0 


(2.9) 


where $ is the implicit function for the node point (i,j). Newton- 
i » j 

Raphson method has been employed here to solve these equations. 
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Based on this method, Q* can be improved by calculating the dif- 

i» J 

ference A Q n which is found from first order approximation of Taylor's 

1 » J 

expansion for $ . Using Taylor's expansion, eq. (2.9) becomes 

i* J 
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Since Q are known from the previous iteration, the set of linear 
i » J 

simultaneous equations (2,12) are solved for the difference A Q, by 

*•» J 

Gaussian Elimination Method. 

The new iterated Q becomes. 


<;- q Lj + aq i.j 


and the procedure is repeated until all A Q. , becomes less than a 

bJ 

prescribed convergence error. 
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A computer-program RSS-FILM has been written to calculate Q for 

i» J 

a given nondimensionaless, gas-film thickness H and a nondimens iona less, 
velocity H. The nondlmenslonal, pressure distribution can be obtained by 
dividing Q by the nondlmenslonal gas film thickness . H. The horsepower 
required to overcome the friction resistance is calculates as 


3 R 0 

w* 0 r 2 3 
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( 2 . 13 ) 


and 


HP = 


n x T x N x 60 
63000 


( 2 . 14 ) 


Also, the load par pad is obtained by the following integration 


W “ p j (P - l)r dOdr 
a J . 

A 


P r 

a o J 


R e r +e. 

, O f> G L 
RdR 


(P - l)d0 


(2.15) 


Figures 3, to 7 show the effects of change In H and H on the pressure 
distributions in a shrouded, Rayleigh- step gas pad. 
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2.4 Approximation for Gas-Film Forces 

The gas- film forces are shown to be dependent upon H, 
following operating and geometrical parameters 
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Shroud Width Parameter 


Length to Width Ratio 


2* 


Table 1 lists the gas-film forces for an annular bearing surface 
containing 20 Shrouded-Ray leigh -Step pads, with geometrical dimensions 
shown in Fig. 1 . By plotting data on a log- log scale in Fig. 8 6t Fig. 9 

it is found that they can be fitted with the following function for a wide 
range of dimensionless thickness, H. 


F 





(2.16) 


Tables 2 list the approximate, gas -film force based on Eq. (2,16) and 
the actual values interpolated from Table 1. The errors introduced by 
using the fitted Equation (2.16) are also listed in these tables. 

The values for C^, C 2 , n^, and n^ for the gas- film forces listed in 
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Table 1 are found to be 


C x « 2.32 lb 

C 2 * 0.76 ^/(in/sec) 

“1- 2.5 

n 2 = 


The above constants, of course, only apply to the operating and 
geometrical parameters shown in Fig.l , For other parameters, a dif- 
ferent set of Cj, C^, n^ and will be required to approximate the 
gas-film forces. 

Letting be the equilibrium nondimens Iona 1 gas film thickness, 
and defining x as the nondimens ionalized displacement from H^, positive for 
a decreasing of H, the gas film force can be alternatively expressed 

as. 


F 





(2.17) 


where H = H - X 

o 


and 


dH 

dT 


dX 

dT 



III. nonlinear axial response 

3.1 Mathematical Modeling 

The main problem in this chapter is to determine the dynamic response 
of the stationary ring to any external force excitation or to any disturbances 
produced by the rotor misalignment and by the rotor surface waviness. 

Knowing the detailed motion of the stationary ring with respect to the 
rotor motion, one may calculate time -dependent gap distribution between the 
surfaces. In general, the dynamic response of the stationary ring is 
measured by the axial translation and by the rotations about two mutually 
perpendicular diameters of the ring. However, for a stationary ring with 
a narrow width and a large diameter, the response in the axial mode is 
very weakly coupled with the oscillation in the angular mode. Furthermore, 
the equation governing the angular oscillation is very nearly the same 
as that governing the axial motion. Thus, it is only necessary to 
concentrate the analysis on the non-linear characteristics of the motion 
in the axial mode. This reduces the problem from a complicated, three- 
degrees-of -freedom dynamical system to a single -degree -of- freedom problem 
for which a more thorough analysis can be afforded. The weak coupling 
between the axial and angular oscillation is demonstrated analytically in 
Reference 1{». 

Figure 10 shows the mathematical modelling of the non-linear 
vibration of the stationary ring in the axial direction. It consists of 
a stationary ring of mass m subjected to a steady load W Q . The back face 
of the ring is flexibly mounted to the frame through a soft spring of 
stiffness and the front face is supported on a very stiff, non-linear 
gas- film whose restoring force is represented by the power relations 
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formulated in Chapter 2, Eq. (2.17). The major problem here is to 
investigate the motion of the stationary ring for one of the two following 
conditions : 

a, a force excitation, q cos ou^t, acting on the stationary ring. 

b, a prescribed rotor disturbance characterized by a Fourier series 
N 

v 

) c cos ncu^t . 

L, n f 

n=l 


3. 2 Equation of Motion 
3.2.1 Force "Excited Notion 

The equation of motion due to a force excitation q cos ajt is 
considered first. Recalling from Eq. (2.16), the force balance on the 
stationary ring gives the following equation. 


d 2 x fl C 2 dh 

"dt 2 H 2 - 5 *H 2 - 5dt 


W q + k s x - q costUft 


( 3 . 1 ) 


where k s x is small comparing with W Q and other terms. 


Introducing the following nondimens ional variables 


X « 

*/6 

T » 

(JJjt 


*f 

U) 35 

* 

B = 

(U 

C 


m6(tu* ) 



where tu is the natural frequency of the system based on the linearized 
equation of Eq. (3.1). Equation (3.1) becomes 


uj?X + CX — — ^ 


(H o « X) 


2.5 


+ B 


f 1 . 

^ (H X) 2#5 H 2 * 5 ^ 
o 


» Q cos T 


(3.2) 


where 


H - H - X 
o 


W 


o 2.5 

o 


d 

dT 


3.2.2 Pis placement -Excited Motion 


N 


cos ntu^t resulted from 


Consider now the disturbance function Y e 

fcl n 

the rotor misalignment, commonly known as the run out, and the non- 
flatness of the rotor-surface. The film thickness H will be perturbed 
by this disturbance, and becomes 


N 

H - H - X + y E cos n T 
o n 

n-1 


(3.3) 


e . 

where E = n/ 6 . 
n 
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The equation of motion in the absence of the excitation q cos 
becomes 


(JD^t 





C 2 dh 
H 2 * 5 dt 


0 


(3.4) 


N 

Substituting Eq. (3.3) into (3.4) and letting X' « X - Y E cos n T , 

n 

one obtains n, “^ 


-2 

cu 


X’ 


+ CX' 


(H - X') 
o 


2.5 


4- B 


I 2 5 

L (H - X') 
o 


H 

o 


1 1 
2.5J 


N 

c- 

/ Q' cos n T 

L., n 

n=l 


(3.5) 


where 



2 

n (d E 


(3.6) 


Equation (3.5) is identical to (3.2) with the exception that X, and 

N 

Q cos ojt are replaced by X' and y Q’^ cos nT respectively. Thus, the 

n=l 

solution for the force-excited oscillations is also applicable to the 

displacement-excited motion provided the proper substitutions are made 

for X' and Q* . 

n 

3.3 Linearized Solution 

If the motion of the stationary ring is such that the resulting 

gap variation, H - H , is only a small fraction of the equilibrium film 

thickness the response can be estimated from the solution of the 

linearized equation about the equilibrium film thickness H . Linearizin- 

o 

Eq. (3.2) for the force-excited motion, one obtains 
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-2 " c 

(ju X H r~c X + X = Q cos T (3.7) 

H ° 
o 

Similarity, linearization of Eq. (3.4) for the displacement-induced 
motion leads to 


—2 X + - 
0) 


2.5 


X + X - e cos T 


(3.8) 


o 

Eqa. (3.7) and (3.8) are clearly the standard, damped vibration equation 
for a single mass, and its solution can be readily written as 


A == 


(1 * «V + S-r t 

H J 
o 


o = tan 


n *^2v 2.5 

. (1 - a) )H 
-1 o 


(3.9) 


(3.10) 


where X = A cos (T - cy) 

3.4 Non-Linear Solution 

Two methods have been employed to obtain the non-linear response 
characterized by the solution to Eqs. (3.2). The first is the method of 
Galerkin.. (Re f . 18) which gives an approximate solution to the nonlinear 
equation. The degree of approximation is governed by the number of terms 
considered in the assumed function in the Galerkin procedure. The second 
method is the direct, step by step, numerical integration using a 
Runge-Kutta procedure. Details of these two methods are given next. 

3.4.1 Method of Galerkin 

The non-linear equation in question, Eq. (3.2), can be represented, 
implicitly as 
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f(X, X, X, T) - 0 

where 

f =* cu 2 X + C * * B ? ^ - Q cos T = 0 (3.11) 

(H - H 

o o 

According to the method of Galerkln, one may assume that the unknown 
response X(T) can be represented approximately by a truncated Fourier 
series , 

N 

X =* ; a cos nT + b sin nT (3.12) 

n n 

n=0 

The substitution of (3.12) into the differential equation gives 
arise to the residue function, 

R(T) = f(X, X, X, T) (3.13) 


The R(T) will not vanish unless X(T) exactly satisfy the differential 

equation. The Galerkin method provides a set of equations by which one 

can solve for the constants a and b for which the residue function will 

n n 

be made extremely small. These conditions are obtained by requiring that 
the integration of the residue function as weighted by each individual 
Fourier component (cos nT or sin nT) be made equal to zero. Stating 
mathematically, one obtains. 



2rr 

R(T) sin nT dT « 0 
o 


( 3.14, 


for n - 0, 1 , 2. , . .N 
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where 

N 

^-2 2 

R(T) =* / *• uo n (a cos nT + b sin nT) 
L. n n 

n»0 


N 

r* 

B + C ' n(-a sin nT + b cos nT) 
i--< n n 

n«0 


N 


H 

L O 


11=0 


(a cos nT + b sin nT) 
n n J 


1 5/ 


2 It 


5/2 


- Q cos 1 


(3.15) 


For N > 1, Integration of Equation (3.14) involves definite integrals 
_ 5 

of - 2 th power of the Fourier series. A gallant attempt was made in 
reducing these integrals in some manageable form, but was unsuccesful. 

Thus, the inclusion of any terms beyond N = 1 was not made. 

For N » 1, 

X = a Q + a x cos T + b x sin T (3.16) 

Further simplification is made by representing X with 


X - A cos (T - o) - A 

o 


(3.17) 


where 


A = - a 
o o 


A - (a. -fb, ) 


2 , 1/2 


O' = 


-i b i 
tan — - 
a. 



21 - 


Substituting (3.17) Into (3.14), one obtains 


,2n 


L 


j | - V0~ A cos T + 


B - CA sin T 


r ' 5 f A l 2 ' 5 

I H + A j ! 1 - „ 2 a cos T 

L o oju H+A J 


- Q(cos a cos T - sin a sin T) J J 


1 

sin T 
cos T 


V dT * 0 


Integrating Eqs. (3.18) and noting the following relations 
,2 tt 


cos T 
sin T 

sin T cos T 


y dT = 0 


.2tt 


2„ 

cos T 

sin 2 T 


dT * TT 


J' 


2rr sin T cos T 


s 2n 


1 - 


V H+A 
o o 


sin T 


cos 


t) 


5/2 


dT = 0 


5/2 


dT = 0 


1 - 

\ H+A 


cos T 


Equation (3.18) reduces to the following algebraic equations 


— 2 H li - Q cos o=0 


- A a) tt 


CA 

- — I. + Q sin 
tt 4 


(3.18) 


(3.19) 


(3.20) 


(3.21) 


(3. 22) 


(3.23) 

( 3 . 24 ) 


2n = 0 


r 3 . :>s> 
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where 


I, = 


_2 ; 

D, J 


S 1 cos T dT 


TJ 

2 

I. = ■— i S. Sin 2 T dT 
4 D- .j 2 

1 G 


2H, 2 ' 5 .2 

*5 = “~dI j S 2 dT 

J. o 


(3.26) 


(3.27) 


(3.28) 




1_ 

A 

H + A 
o o 



/ 

. 1 


H + A 


cos 


_ v 

T ) 

/ 


2.5 


(3.29) 


D 


1 


(H 

o 


+ A ) 
o 


5/2 


(3.30) 


Eliminating o’ from (3.23) and (3.24), one obtains a system of two non- 
linear equations to be solved for A and A . These two equations are 


/ BI l -ox 2 /CAI. 2 _ 

F(A -V = ("T - A ■ ) + (it) - Q - 0 

G(A,A ) - I_ - 2tt = 0 
o 5 


(3.31) 

(3.32) 


Using Newton-Raphson procedure 
be expressed in terms of F, G, 
last inter ate s of A and A . 



o 


3F 

-F 

dA 

o 

-G 

dG 


dA o 


, the successive corrections A A. A A can 

o 

dF dF dG dC 

SA* sT’ 3A’ and aT’ evaluated at the 


(3.33) 
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A A 


where A = 


a f 

-F 

dA 

. I 

ac 

-G 

dA 


9F 

dF 

dA 

dA 


o 

dG 

dG 

dA 

dA o 


( 3 . 34 ) 


( 3 . 35 ) 


5F _ / BI 1 
dA ~ \ TT 




ai, 

ar) 


(3.36) 


JE- . 2 £i 

dA V 

o 


\ TT 


ai, 2 2 ai. 
A ) n aT + 2 ^2 *4 aT 


( 3 . 37 ) 


dA 5A 


(3.38) 


dG 

dA 

o 



(3.39) 


The differentials of nonlinear integrals can be derived from equations 
(3,36) - (3.28). They are 




TT 

9 


cos T S dT 
o 


} 




cos T S^dT 


} 


( 3 . <40 i 



24 


Equation (3.40) cont’d, 


rr 


51 4 _ 1 r 5 

dA D 2 l 


sin^ T cos T S^dT^ 


dl 


5T “ S i 5 1 sin2 T s 3 dT ’ 

O l o 


n 


ai, h 2,5 r .2 

_5 . { 5 j^cos T S/Tj 


n 


dA~ 


2 5 

H ^* ** f .2 


V i 5 j o s *"} 


where = (H + A ) 

z o o 


3.5 


and 


/ 

f 1 

\ 


H + A 
o o 


cos T 


3.5 t +} , 

) c 


3.5 


1*1 + 


H + A 
o o 


cos T 


Being provided with values of A, and A^, a can be solved from Equations 


(3.23) and (3.24) as 


a' ~ 


tan 


-1 


CAI, 


BI^ - tt A 


O . 4 i > 
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A computer program has been written to solve for A, A and o. from which. 

o 

X and X . can be determined by 
max min 3 


X 

max 

X min 



(3.42! 


The integrals 1^,, 1^, 1^ and their derivatives with respect to A and 
A^ are evaluated numerically. Tables have been prepared for various 
values of A/(H q + A q ) varying in the range, [0,1] at intervals of 0.01. 

3,4.2 Direct Integration 

The step-by-step numerical integration of Eq. (3.2) is achieved by 
splitting it into two first order equations: 


X = Y 




B 


cu 2 l 'H 2 ‘ 5 (H - X) 2 ' 5 
o o 


B + CY , _ ^ 

+ Q cos Tj 


( Ref. 16 ) 

The popular Runge-Kutta method has been employed in obtaining the 

solution for X and Y with a given set of initial conditions X , Y . The 

o’ o 

response is represented by the trajectories in the phase space plot (X,Y). 
The response is considered to have reached a steady state if the trajectory 
approaches a limit cycle, which could be a single or multiple- looped cycle, 
A library subroutine at the Vogelback computer center has been used for 
this numerical integration. Fortran listings for the Computer program, 
RSGALN, which calculates A, A q , and by the Galerkin method, and the 
computer program, RSRKIT, which calculates the trajectories in the phase 
space, are included in the Appendices and 
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3.5 Results of Non-Linear Response 

The results of the nonlinear response are presented in two parts. 

The first part is obtained from the method of Galerkin one -term approxi- 
mation, and the second part from the Runge-Kutta direct integration. 

Since the gas-film is an unsymmetrical spring, i,e, , the relation between 
the displacement and restoring force is not symmetric with respect to 
the equilibrium position, the amplitude of response during an upstroke 
is different from that during a downstroke. Duriug an upstroke, the 
gas film stiffness is softer, and the amplitude is greater than that during 
a downstroke when the gas-film is considerably sttffer. For this reason, the 
response in the upstroke and downstroke are plotted separately against 
the non- dimensional excitation frequency in Figs, 11 to 19 • 

3.5.1 Results by Method of Galerkin 

Referring to the equation of motion, Eq„(3.2',it is seen that the 
parameters affecting the dynamic response are: 

B = stiffness parameter — 

C = damping parameter = 

= static-film parameters = 

Q = forcing intensity parameter - 

tu = frequency parameter - 


* 2 

m6(cu ) 


C 2“f 


mo; 


7 6 


* 2 
m6(o) ) 




f/ * 

09 


The general characteristics of the upper and lower amplitudes as 
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a function of excitation frequency are illustrated in the sketch below. 



The nonlinear response based on one -term Galerkin method is shown as 
solid lines except with a small portion of the unstable oscillations , which 
are shown as dotted lines. The dashed lines show the response predicted 
from the lin arized solution. In the region AB, the excitation f requeue- 
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is smaller than the natural frequency, u) , based on the linear theory* 
and the non-linear solution predicts a smaller lower amplitude but a 
greater upper amplitude comparing to the linear response. As increases 
to the region BC, the non-linear theory yields three possible solutions, 
one along the path BC, one along the path B’ C, and another one along the 
path B'C'. The solution along BC is in-phase with the forcing function 
and is the most stable mode of response; the solution along B’C is un- 
stable and only exists mathematically; and the solution along B’C 1 is 
out-of-phase with the excitation and is less stable than the solution along 
BC. For excitation frequency beyond the region BC, the characteristics 

of the non-linear response are similar to that of the linear response in 

* 

the region where uj^ > u) . In this region, the pad is insensitive to 
the excitation and would not track any disturbance introduced by rotor 
runout or waviness. 

The non- symmetrical nonlinear gas-film produces a response charac- 
teristics which resemble more to the response due to a symmetrical, 
soft, nonlinear spring, for which the resonance occurs at a frequency 
considerably lower than the natural frequency based on the linear theory. 
This correlation is really not surprising, since the mean position of the 
oscillation shifts to the region of softer stiffness, and the nonlinear 
oscillations are dominated by the softer part of the gas-film stiffness. 

Fig. ll shows typical response curves for the following dimensional 
parameters : 

h = 0.0005 inches 

o 


6 = 0.001 


inches 
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m = 0.2 slug or 6.44 Ibm 

C x * 2.32 lb 

C 2 - 0.76 lb / (in /sec) 

q =1,2, and 5 lb 

The corresponding ,non-d intension a l parameters are listed in Fig, 11 . 

The inward bending of the resonant peak in the region id < 1 is clearly 
visible in alL three cases. Fig. 12 shows the effect of increasing the 
mass of the pad from 0.2 slug to 1.0 slug. The increase in mass does 
not alter the parameters, B, fi^, and Q, since (cu ) is inversely pro- 
portional to m. The only parameter affected by changing of m is the non- 
dimensional damping parameter. A five fold increase in mass is equivalent 
to a ^ tiroes reduction in the effective damping factor C^. The more 
peaky response near the resonance is clearly visible in Fig. 12 when the 
mass is increased by five fold. 

Figs. 13 and 14 shows the effect of increasing the equilibrium film 
thickness from 0.5 to 0.75. The natural frequency is reduced sharply by 
the increase in the film thickness, and the level of response is also 
much greater with a thick film than with a thin film for the same forcing 
function. 

To investigate the effect of damping the value of C 2 has been doubled 
and halved from the case shown in Fig. 11 . The curves in Fig. 15 show 
that when the damping is doubled, the response near the resonance is 
considerably suppressed. The opposite effect is introduced if the damping 
is halved as shown in Fig. 16 . 
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3 , 5*2 Results by Direct Integration 

Both the upper and lower amplitudes obtained by using the step-by- 
step, Runge-Kutta, direct integration are plotted against the excitation 
frequencies in Figs* 17 and 18 . A case of heavy mass, small equilibrium 

film thickness, and large excitation force has been selected to illustrate 
the nonlinear effects. The linear response curve and the approximate 
nonlinear response by Galerkin method are also plotted as dashed and 
dotted lines for comparison. It is seen that the agreement between the 
Runge-Kutta results and the Galerkin results is good near o> = 1. This 
clearly shows that even with one term approximation the Galerkin method 
yields a reasonably accurate prediction for the synchronous response. 

For u) < 1 , the Runge-Kutta results show a series of superharmonic reso- 
nances at ui approximately equal to 1/2, 1/3, 1/4, etc. The magnitude of 
these super harmonic amplitudes is, of course, governed by the damping 
factor. 

Fig. 19 shows the trajectory in the phase-space plot for condition 
near the second superharmonic resonance. The final limit cycle forms a 
two-loop orbit showing typical characteristics of a superharmonic response. 
Other trajectories at the third and fourth superharmonic resonances are 
also shown in Figs. 20 to 22 . A subharmonic resonance is also found 

for u> — 2.0, but the amplitude is small and harmless. The Characteristics of 
the phase space trajectories near uj = 1 are plotted in Figs. 23 to 25. 



IV - STABILITY OF AN INFINITELY- WIDE RAYLEIGH- STEP FAD 
4.1 Statement of the Problem 

It is well known in hydrodynamic lubrication that a dynamic system 
involving any fluid- film supports may, under certain conditions, encounter 
detrimental self-excited oscillations commonly known as dynamical instability 
of fluid- film bearings. The gas-film bearings are particularly susceptible 
to this type of instability. The fractional frequency whirl of a 
shaft supported on gas-bearings and the pneumatic hamner in externally 
pressurized gas-bearings are two of the prominent examples of the fluid- 
film instability. For journal bearings, the gas-film instability usually 
occurs if either the rotating frequency or the supported mass becomes 
large. There have been considerable data available to predict the threshold 
speed or critical mass of the journal bearing. However, 

for gas- lubricated thrust pads, the problem of instability is relatively 
unexplored. Since present trends in gas-bearing are always toward higher 
and higher speeds, it is important to determine whether there exists any 
stability threshold associated with a gas- lubricated^ thrust pad. 

This chapter is devoted to the stability analysis of a thrust pad 
with a Rayleigh -step. The geometry of such a thrust pad is shown in 
Fig. 26 . In order not to impose excessive burdens on the analysis, the 
assumption of an infinitely wide pad has been made. Moreover, the motion 
of the pad Is assumed to be restricted in the transverse direction only. 

With these two assumptions, the problem is reduced to the stability of a 
single-degree-of- freedom dynamical system with restoring pressures 
governed by a partial differential equation in space and time. 
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4,2 Governing Equations 

The one -dimensional, time -dependent pressure distribution is governed 
by the Reynolds equation, 


iL 

dx 


( ph3 "I) = 6m “ 


$ishl + 12u 
ax * at 


(4,1) 


where h - 6 + + e(t) for 0 < x < B^ 

h 33 + e (t) for < x < B 

e(t) 1,3 the upward motion of the pad. 

The boundary conditions for Eq. (4.1) are 
p = p at X * 0 and x * B 

rl 

The equation governing the pad motion is. 



where m is the mass per unit length of the pad, B the length of the pad, 
and W q the static load imposed on the pad. The dynamic system represented 
by the coupled equations, (4.1) and (4.3) are to be investigated for the 
stability for an equilibrium position. 

4.3 Method of Solution 

The time -dependent, nonlinear, pressure equation, Eq. (4.1), is 
difficult to solve analytically, A numerical approach has been used here 

in solving a set of discretized, time -dependent, pressures along the 
coordinate X. 

The discretization of pressure is achieved by considering the flow 
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balance in an elemental volume within the gas film as shown, in Fig, 26 for 
the i-th element. The flow rates in and out the element are respectively 
q j and and 


<1 


«2 



(. J*L SE + U 1 + U 2 

V 12p dx 2 



where 


u. 


u and = 0 


(4.4) 


Considering the flow balance, 

^1 ~ ^2 ** ra ^ e “ass stored within the elemental volume 
It follows. 



C* 




+ h 


i+% 



) 


n 

1 

J 


(4.5) 


Using the isothermal relation. 


~ R T = constant 
p e 


Introducing the non-dimensional parameters. 



X -» X /B 


(4.6) 


( 4 . / ) 



Equation (4.7) cont'd 


H = k - h 


6 h 2 - h x 


A « 


6uBU 

c2 

P a 6 


12(j,vB' 

(j SI * — 

f .2 

P o 
a 


T « vt 

V excitation frequency 


and using the following finite approximations for ; p H 

\SX/ » r i-Y H i-j 


,M j _ !l P i ' P i-1 

*i - Vl A 


i-% 


2 <P i + P i-1> 


‘i-% “ 2 (H i + H i-l) 


the discretized pressure equation becomes 


(4.8) 


A-_ / P 2 . 2^ 


4A X 4 l p w - Vj H s,i - 4^1^ ( P i 2 - P 4 .J) H Sii _ a 


‘ A (P i+l + P i )(H 1+ l + «!> + A (Pi + Pj^HHi 4 H.^) 


= a y [p t (H t + H^)] A X 1 . 1 4 y ;> i( H i+1 4 V 


a x, > 

ij 


where 


H . » i H * ” 


(4 
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where A « a .P, - + d . 

i o,i 1+1, o o,i 


B = b P + e 
i o, i i,o o,l 


C ■* e P - + f 

l o,i i-l,o o,i 


H 


S.JL 


0,1 4 A X, 


, d * - A(H , + H ) 

o,i l+l,o i,o 


H 


o.,i 4A X 


3fl" 3- 


1-1 


, f . = A(H + H. . ) 

o,i ' i,o i-l,o 


b o : ,i - - a o,l - C o,i > e o,i “ d 0 ,l + f o,i 


(4.17) 


and 


P. = P - 1 
l,o n,o 


(4.18) 


Since the algebraic equations (4.16) are nonlinear, Taylor's expansion 
is used to reach the following simultaneous equations for A P 


A *1 " “ V P i,o> 


(4.19) 


where A * = (A + a .P.._ )A P, Lt 

1 1 o , i i+l,o i+l,o 


+ + b o J A 

i 0,1 1,0 1,0 


+ (C, + C .P. ) A P. - 

1 o,i i,o i-l,o 


Equations(4 , 19) are inverted directly for successive A P until the 

i,° 


convergence is reached 
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In equation (4,9), both the discretized pressure P^ and are dependent 
on time. For transient studies, they have to be solved simultaneously with 
the dynamic equation of motion, Eq. (4.3). However, for small oscillations 
and stability analysis, the variation of H and P with time can be con- 
sidered as small perturbations about the equilibrium solution, H q and P^. 
These small perturbed quantities can be expressed as 

H^T.X) - H (X) + e(T) (4.11) 

P (T,X) = P (X) - e(T) P (X) (4,12) 

i i,o i,c 

where e * e /6 

l«il « 1 

The minus sign in (4.12) is due to the fact that an increase in 
thickness leads to a decrease of pressure of the gas film. Substituting 
(4.11) and (4.12) into Equation (4.9), the equilibrium equation and the first 
order^perturbed equations can be obtained. They are 

4A X ± V'i+l,o *i,oJ s,i 4A V i,o P i~l,oJ s,i-l 


- A (P.,* + p, ) (H . + H ) + A (P, + P. , )(H. + H. 

i+l,o i,o i+l,o i,o i,o i-l,o i,o i-l,o 


SB 0 


(4.13) 


and 


(P l+l,o a s,i + a n,i )P l+l,c + 


( p , b . + b )F 
i,o s,i n, i i,c 


(P, . c . + C .) - d, 
i-l,o s, l n,v i 


1 de T 

e AT [. 


P. 

1,0 


f. 

l 


+ e.P i 
i i,c_j 


1 4 . 14 i 



- 37 - 


where 


H 


8 jlL 


s,i 8 A X J 


H 


b j1zL 


s,i 8 A X 


i -1 


\ * - 7 (H + 1 + ) 

n , i 4 i+1 , o i , o 


c » 7 (H. + H 4 . ) 

n,i 4 i,o i~l,o 


b * •• a . - c . , b , = a . + c 

s,i s,i s,i 9 n,i n,i n,i 


d i - XT ( P m,o 2 - P i,o 2 )( H i + l,c 2 + H i + l,o H i,o + H i,o 2 ) 

- ( P i,o 2 - P i-l,o 2 X«i,o 2 + H i,oVl,o + H i-l,o 2 ) 


2 (P i+l,o ’ P i,o> 


». - 7 ! (H. - + H ) A x. + (H + H ) A X j 

i 4 l i+l,o i,o 1 i,o i-l,o 1-1J 


£ i - - f X i-1 + A X i> 


(4.15) 


Equations (4.13) are nonlinear simultaneous algebraic equations for 
the pressure distribution P^ q . They are solved numerically by Newton - 
Raphson method. The boundary conditions are P^( l*n ) *= 1 , 

Rewrite equation (4,13 ) as 


> 

i i , o 


A P _ + B , p , + c p 4 i = 0 

i i+l,o i i,o x i-l,o 


for i «* 2, . . . ,n - 1. 


(4.16) 
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After the equilibrium pressure distribution is found. Equations 

(4.14) are solved for the first order pressure distribution, P . 

i,c 

In Equations (4.14), both ^ and e(T) are complex quantities, and they 
are assumed to be 


P * P + IP 
i,c i,R jr i,I 


and 


e ~ e e 

o 


j T 


(4.20) 


(4.21) 


The boundary conditions are P « 0 at entrance and exit, which gives 

l , c 


P M eP n,R sP l,r P n,I e0 


(4.22) 


Substituting Eqs. (4.20) and (4.21) into Eq. (4.14), one obtains a 
system of 2n ~ 4 simultanous, linear equations. 


A i+l,i X i+l + A i,i X i + A i-l,i X i-l + A i+ n -2,i X 14n-2 * D i 

for i“ 1, ••••••• n _2 


(4.23) 


"1-^,1 l-nW + *1+1,1 *1+1 + *i,i X i + *i.i,i X i- X - »i 

for 1 • n-1, f 2n-4 

In the above, A and X r— * ** i . . 

* «• ana a for x - i, n-2 are defined by, 

j - 1+1 


A i+1, i ” P j+l,o a s,j + a n,j 

’ X t+1 “ P *«,R 

A. » P. b . + b , 

i * 1 j.® n,j 

» x. « P 

1 i+l,R 

i“l» i s,j n, j 

, X, , - P 

1-1 1,R 

A i+n-2,i " e j 

• X i+n-2 =P l+l 


(4.24) 


D. = d 

1 J 
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and for i = n - l,...,2n - 4, they are defined by, 
j * i-n+3 


i+l,i P j+l,o a s,j + a n,j 

' X i+1 

« p 

i-n+4,I 

, , = P. b . + b . , 

i* 1 s,j n,j 

■ X l = 

P i-n+3,l 

u . - P. , C + C . . 

i-l,i j-l,o s,j n, j 

• X i-1 

" P i-irf2,I 

= ' e j 

’ X i-n+2 “ P i-tr+3 

* f P 
1 j J »° 




The boundary conditions become 


(4.25) 


X 

o 


1,R 


0 


X 

n- 


1 



for equations with i = 1, 
n - 2, respectively 


and 


X n-2 - P 1,X * 0 


x. = p = o 

2n-3 n,l 


for equations with i = n - 1, 
2n - 4 respectively 


(4.26) 


(4.27) 


These simultaneous equations (4.23) can now be solved for fX] vector 

by direct matrix inversion. A computer program has been written to solve 

for P, and P for different real values of v, and the Fortran listing 

of this program is included in Appendix D. 

Once the real and imaginary part of the gas film pressure are determined, 

the integration of P. and P. gives respectively the in-phase and 

1 , K 1) i 
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,1 


out-of -phase bearing forces. The in-phase force { P dX, can be interpreted 

J o R 1 

as the stiffness of the film, whereas the out-of-phase, j P dX, represents 

J o 

the damping factor of the film. It should be emphasized that both the 
in-phase and out-of-phase perturbed pressure are dependent upon a, or 
the excitation frequency v of the mass. The frequency-dependent charac- 
teristics of the gas- film reactions is a direct consequence of the inclusion 
of the term ~ . The use of these frequency-dependent bearing forces in 
determining the dynamic stability of the gas-film and pad system is des- 
cribed in the next section. 


4.4 Stability Criterion 

The stability of the gas-film and pad is governed by the equation of 
motion, Eq. (4.3), which, in its non-dimensional form, appears as 

j2 p B 1 W 

~ 2 “ ~ S ~2 (P ' 1 > dX ' (4.28) 

dT m6v o m6v 


Recalling the pressure is the summation of the equilibrium pressure, P , 
and the dynamic pressure, - eP^, one obtains 


p B .1 W p B 1 

-Z-J J (P - l)dx = . _a | p dx 

m 6v o m 6v m &v o C 


(4.29) 


It follows that 

d 2 1 

— f + e I P dX = 0 (4.30i 

dT J o C 

Mathematically, Eq. (4.30) represents a free oscillation problem which 
contains stiffness and damping factors depending upon the frequency of 
the oscillations. A direct approach in determining of the stability of this 


; 'm6y 

Vp B 
a 
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single-degree of freedom problem is to look for the eigenvalue of this 

system. If the real part of the eigenvalue is negative , the system is 

stable; otherwise it is unstable. .If the eigenvalue is purely imaginary, 

then the system is at its threshold of instability. 

Alternatively, one can also determine the stability thresholds by 

assuming that the eigenvalue is a purely imaginary number and inquire what 

mfiv 2 

would be the mass parameter, — g~, for a purely imaginary eigenvalue. 

Let the eigenvalue be represented by K + jv# and for a pure Imaginary 
eigenvalue, \ = 0, It follows that 


jvt jT 

e - e e J * e e J 
o o 


(4.31) 


Substituting Eq. (4.27) into Eq. (4.26), and separating the real from the 
imaginary part, one obtains 


„1 


J PjCvJdX « 0 


(4.32) 


j 1 p R(v)(lx - 2^ - 0 
o a 


(4.33) 


where Pj(v) and P R (v) are solutions of Equations (4.23) for a given value 

of v. The pure imaginary eigenvalue v may be determined by evaluating the 
1 

integral, j P (v)dX , for various values of v until the integral changes 

its sign. The exact eigen frequency v may be calculated by a linear 
interpolation. Once the eigen frequency is found, the critical mass at the 
threshold of instability can be determined by Eq. (4.28), and 

P a B .1 

m cr = 75 J V v)dX 

6v o 
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Equation (4.29) predicts a quantitative value of the critical mass, 
but does not furnish any information on ’which side the stable region lies. 

To determine the region of stability, one may use the criterion developed 
by Malanowski and Pan (ref, 17). Their criterion can be stated in the 
following manner. 

Stability Criterion - The system consisting of a thrust bearing of mass, m, 
with Rayleigh Shrouded-Step Seal is in a state of self-sustained oscillation 
at frequency v when and only when the out-of-phase component of the 

° r 1 

bearing reaction, P (v)dX, vanishes, and when the mass, m has the 

j i. o 

o 

critical magnitude to be in resonance with the in-phase bearing reaction. 
Equation (4.28). The system will become unstable if the mass exceeds the 
critical value provided the out-of-phase bearing reaction increases with 
the frequency in the algebraic sense in the neighborhood of the critical 
frequency V 0 ? and conversely. 

4.5 Results 

The steady state pressure distribution, P , has been calculated 
numerically for the following parameters: 

®1/ B = 0.5, and 0.75 

H » -■ = 0.5 and 0.75 

0 

A -= " ~ 8.4, 42, and 100 

The resulting pressure curves are plotted in Fig. 27 , and they are in 

excellent agreement with the analytical solution provided by Kochi ( R e f. 6) 
Chis comparison confirms the accuracy of the present numerical solution 
of the steady state pressure distribution,? . 
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For each steady-state pressure distribution, the dynanic pressure 
distributions, P R and P^. are calculated for a series value of a. Figs. 28 
to 35 show a typical series of the dynamic pressure profiles for B^/ B * 0,5 

and 0.75 H * 0.5, and A = 42. For small excitation frequencies, the 
real part of the dynamic pressure is dominated by the bearing parameter A, 
and the profile is similar to the static pressure distribution, and is 
relatively independent of u, As a increases, the pressure distribution, 

Fp > becomes slightly wavy at both edges. The waviness penetrates deeper as 
a further increases. As a approaches infinity, the effect of A disappears, 
and P R approaches the asymptotic solution, 

P 

(P ) = — 

R Q - ® R o (4.35) 

The imaginary part of the dynamic pressure takes a wavy pattern even for 
small values of <j. For a approaches infinity, the values for P vanishes 
throughout the entire region. 

The in-phase and out-of-phase forces are plotted as a function of 
a in Fig. 36 , and also listed in Tables 3 to 5. It is seen that for small 

or moderate values of A, the out-of -phase force never becomes negative. 

This indicates that for nearly incompressible cases, there exists no stability 
threshold, and all equilibrium solutions are stable. As A becomes ex- 
tremely large, the out-of-phase force does become negative at a fairly 
high value of a. This crossing-over of zero line indicates that for a 
highly compressible film, there does exist a stability threshold, and the 
gas film will exhibit a self-excited oscillation at a fairly high frequency. 
Moreover, the criterion in Reference 16 suggests that the stable region 
lies in the area where the mass of the pad is greater than the critics; 
mass calculated according to Eq. (4.34). 



V, SUMMARY OF RESULTS 

1. The gas-film rescoring forces in a Rayleigh- Shrouded- Step thrust 
pad can be determined numerically by solving the discretized time-dependent 
Reynolds equation with irregular grid-spacings to account for any abrupt 
changes of pressure at the step and at the exit edge* 

2. For conditions of high ambient pressure, for which the term of 

h can be neglected in comparison with other terms at the right side of 
ot 

the Reynolds equation, the gas- film force is found to be approximately 
inversely proportional to nth power of the film thickness and directly 
proportional to the squeeze-film velocity. The exact value of n is a 
function of the step geometry. In general, n lie between 2 and 3. 

3. The axial, non-linear response of the Rayleigh- Shrouded-Step 
pad to a sinusoidal, axial forcing function or a sinusoidal disturbance 
due to the rotor misalignment or surface waviness can be determined by 
one of the following two methods: 

a. By assuming the response to be a truncated Fourier series 
in multiples of the excitation frequency, the Ritz-Gaierkin 
procedure can be employed to predict the non-linear behavior 
of the pad motion. 

b. By integrating directly the equation of motion of the thrust 
pad using a step-by-step, numerical routine, the Runge-Kutta 
procedure. 

4. Results obtained by using the Ritz-Galerkin method with the 
first harmonic terms show considerably departures from the linear response 
curve as the frequency approaches the resonance based on the linear theory. 
The asymmetric spring characteristics of the gas film result into a non- 
linear response similar to that caused by a syranetric, soft, non-lineai 
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spring. The resonating peak bends Inward and occurs at a frequency leas 
than the resonating frequency based on the linear theory. The peak can 
be suppressed by decreasing the mass. Increasing the damping, and Increase 
the stiffness, 

5. Results obtained by the step-by-step direct Integration confirms 
the approximate solution in the vicinity of the resonance. The direct 
integration also predicts a number of additional peaks at frequencies less 
than the resonating frequency known as the superharmonic resonance. 

6. The gas-film instability of an inf lnitely-wide Rayleigh step 
thrust pad canbe determined by solving the complete, time dependent, 
Reynolds equation coupled with the equation of the motion of the pad. 
Results show that for bearing numbers, A, up to 50, the Rayleigh step 
geometry is very stable, and no stability threshold has been discovered. 
For ultra high values of A i 100, a stability threshold Is shown to exist, 
and the stability can be achieved by making the mass heavier than the 


critical mass 
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TABLE 1 


GAS FILM FORCES (lb f ) 

d - 6.46 Ins d « 5.96 ins 
o i 

6 » 0.001 in (step) 
u> * 277 rev. /sec. 

P * P • 315 psla.: Ambient Pressure 

a o 


^sjmin 

dh/dbv^ 

2.0 

1.5 

1.0 

0.75 

0.50 

0.30 

0.20 

-1 in/sec. 

0.50635 

1.10134 

3.2606 

6.9965 

19.949 

65.138 

139.729 

-0.5 

0.45104 

0.97913 

2.8878 

6.1765 

17.499 

56.041 

115.531 

-0.25 

0.42338 

0.91801 

2.7015 

5.7667 

16.275 

51.529 

103.734 

+0.25 

0.36806 

0.79580 

2.3290 

4.9477 

13.833 

42.576 

80.750 

+0.5 

0.34040 

0.73470 

2.1428 

4.5384 

12.615 

38.135 

69.559 

+ 1 

0.28509 

0.61250 

1.7704 

3.7204 

10. 182 

29.325 

47.778 

0 

0.3957 

0.8569 

2.515 

5.357 

15.06 

47.124 

92.8 

F - F 
gl g2^ 

h 2 - 

0.1106 

0.2444 

0.7450 

1.6380 

4.8836 

17.907 

45.975 


6 x 10 


_9 lb^- sec 

’ ~ 77 ~ 


\ * 6.30076 
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TABLE 2 


LOAD (lb f ) CALCULATED AND ERROR OCCURRED (%) 


VWn . 

h in/sec\ 

2.0 

1.5 

1.0 

0.75 

.50 

.30 

0.20 

-1 

0.544 

1.116 

3.08 

6.33 

17.403 

62.7 

171.2 


+ 1.5% 

+ 1.45% 

- 5.53% 

- 9.5% 

- 12.5% 

- 3.7% 

+ 22.57. 

-0.5 

0.477 

0.978 

2.70 

5.55 

15.28 

54.95 

150.1 


+ 5.78% 

- 0.1% 

- 6.5% 

- 10.1% 

- 6.29% + 6.65% 

+ 30% 

-0.25 

0.4435 

0.909 

2.51 

5.16 

14.205 

51.075 

139.05 


+ 4.78 % 

- 0.98% 

- 7.07% 

- 10.5% 

- 12.3% 

£ 

00 

• 

0 

1 

+ 33% 

+0.25 

0.3765 

0.771 

2.13 

4.38 

12.055 

43.325 

118.95 


+ 2.34% 

- 3.01% 

- 8.2% 

- 11.5% 

- 12.9% + 1.76% 

47.4% 

0.5 

0.344 

0.703 

1.94 

3.99 

10.98 

39.45 

107.9 


+ 1.18% 

- 4.22% 

- 9.35% 

- 11.8% 

- 12.9% + 3.46% 

55.2% 

1 

0.276 

0.564 

1.56 

3,21 

8.83 

31.7 

86.8 


- 3.15% 

- 7.85% 

- 11.85% 

- 13.7% 

- 15% 

+ 8.2% 

+ 82% 
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TABLE 3 


DYNAMIC BEARING REACTION FOR 
B^/B « .75 H 2 - 0.5 
A * 8.4 


r r 

a I FjdX j P R dX 

20 .13373 1.1836 
50 .21003 1.2694 
100 .27600 1.4097 
200 .24248 1.5930 
300 .18121 1.6590 
450 .12761 1.6887 


1 

.005975 

10 

.057536 

50 

.11984 

100 

.020667 

200 

.031105 

300 

.26865 

500 

.37916 

2 x 10 3 

.079279 

10 4 

.045201 


ax 42 

1.5018 

1.5144 

1.6844 

1.7192 

1.5013 

1.5269 

1.8392 

2.0095 

2.0296 


A = 100 


10 

.05897 

1.0295 

50 

.28043 

1.1019 

100 

.47831 

1.3014 

200 

.47977 

1.7884 

300 

.13904 

1.9829 

500 

-.12737 

1.5789 
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TABLE 4 


DYNAMIC BEARING REACTION FOR 
B^B * .7 5 H 2 - 0.75 


A « 8.4 



l* 1 

r i 

a 

i PjdX 

o 

J V* 

o 

10 

.15200 

.55179 

50 

. 16382 

.78737 

100 

.15032 

.84853 

200 

.14273 

.92439 

300 

. 12283 

. 96835 

500 

.08910 

1.0045 


A « 42 


1 

.009264 

.61777 

10 

.090506 

.63251 

50 

.27479 

.86727 

100 

.17722 

1.0606 

200 

.052296 

1.0238 

300 

.10195 

.99666 

500 

.15700 

1.0797 


A «= 100 


1 

.004522 

.46999 

10 

.045139 

.47216 

50 

.21611 

.52320 

100 

.37708 

.66548 

200 

.42895 

1.0325 

300 

.22267 

1.2315 

500 

- . 036699 

1.0725 
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TABLE 5 


DYNAMIC BEARING REACTION FOR 
B 1 /B « 0.5 H 2 * 0.5 

A = 8.4 



A 

p i 

a 

I V* 

J V 


o 

o 

1 

.04398 

.9242 

10 

.40067 

1.0473 

50 

.65404 

1.9523 

100 

.34743 

2.2234 

200 

.17872 

2.2381 

300 

.14173 

2.2543 

500 

.10469 

2.2785 

800 

.07845 

2.2960 


A * 42 


1 

.00519 

2.3313 

10 

.05146 

2.3349 

50 

.22100 

2.3947 

100 

.38003 

2.4791 

200 

.66419 

2.7731 

300 

.65252 

3.1907 

500 

.22427 

3.3674 


A * 100 


l 

.00727 

1.9138 

10 

.07261 

1.9170 

50 

.34848 

1.9943 

100 

.61342 

2.2089 

200 

.73899 

2.7584 

300 

.49508 

3.0561 

500 

.29231 

2.8927 
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r 0 =3.23 in 
ri=2.98 in 

8=0.001 in (step) 
u)=l740 rad/sec 

i£WV a08 

p fl =p o =3l5psia 

A= 8.30076 

0 G = 109° 

0 L =5.32° 



View A-A View B-B 


Fig.l The geometry of a shrouded, Rayleigtrstep 
thrust bearing 
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Fig-2 Flow bolance around atypical grid point 
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Fig. 5 Pressure distribution for h- - 1. i ps. I. 
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Fig 6 Pressure distribution for h»l.ips>Wi 
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Fig. 7 CONTOUR . MAP FOR PRESSURE DISTRIBUTIONS 
for li =-l in/sec, KsiifO.5 
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Fig 8 Variation of static, gas film force (H=0) 
with the normalized film thickness 




Fig.9 Variation with the normalized film thicknes 
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0 . 5 ’ 

I 

Fig.ll 



Fig.12 Nonlinear response for Htf=0.5, m s lslug 







* 1 *) 


w =436 rad/sec 


Q=0.I27 


0=0.0634 


V Q = 003I7 




I 0*0127 
f/ Linearized 
Prediction 


Fig.14 Nonlinear response for Ho=0.75, m=l slug 
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Fig. I 5 Nonlinear response for 


Ho- 0.5, m« 0.2 slug.C -1.52 H^n^ec) 
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u/»l984 rad/sec 



Fig. 16 Nonlinear response for H»»05 
m -0. 2 5 slug C 2 = 038 Ibf/ips 




Q* 0.076 


Ho 0.5 
m= I slug 

+ Runge-Kutto 
— Ga I ertfin Method 

Linearized 
Pre d ic ti on 

887 rad/sec 



Q=a076 

H*05 

m=lslug 

+ Runge Kutta 

- Golerkin Method 

— Linearized 
Prediction 

•887 rad/sec 


Fig 18 Comparison of the approximate and exact downward amplitude 

of response 
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d harmonics 
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Fig.2l Phase plot 
a>= 0.248 


orth harmonics 



Fig.22 Phase plot with 2nd order 
subharmonics 

w = 2.096 



c uj t - 0. 76 6 / ) 
cu — o. 7 3 u 


Fig.23 Phase plot with limiting cycle of 
large amplitude 
£> =0-732 
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Fig.24 Phase plot with limiting cycle 
of small amplitude 
5=0.745 







■h*i ' 'i+ 
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Fig.26 Flow balance and geometry of a infinitely 
wide Rayleigh step pad 
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Fig.27 Pressure profile for B,/B 3 0.5, H,=0.5 


















Fig. 36 Variation of dynamic bearing forces with the excitation 

frequency, <J for 13/13=0.75, A=42, H 2 =0.5 ? 
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PROGRAM RSE AL ( INPUT » OUT PUT • TAPE 5 = INPUT , T APE6*UUTPUT ) 

C RAYTH 

D 1 MENS I ON D X ( 1 7 i *DTH < 33 ) »HM1 NA ( 10 ) * RPSA ( 10) *PAA ( 1Q)*VISA(1Q) *TH< 
1) » R ( 1 7 ) > F RC ( 17 ) » H ( 1 7 * 3 3 ) *A1(17*33)*A2U7*33)*A3(17*33)*A4( 17*33) 
2 A5<17*33) * A6 (17*33) »A8( 17*33)*d(17)*C(17)*Atl7*33)*F( 17) »QSMA( 1 
3l7)*E(17,17*33),G(17*33)»OQli7»33)*P(17*33)*U(17.33) »OQQ{ 33 ) * 

4PP (17) » HOU T ( 10 ) » H DO T (10 ) 

1 FORMAT (72H 

1 ) 

2 FORMAT (1615) 

3 FORMAT ( 8F10. 7) 

+0-MAT (+5*5X,?E10.4/ ( 8E10C4DD 

5 FORM£ T ( 6X * AHlAST » 6X * 4HNPAD * 9X * 1HM * 9X* 1HN*8X# 2H1FU6X, 2HJH*7) 

1 3HIHH.7X* 3 H wi H H * /) 

6 FORMAT ( 81 10 */ ) 

7 FORMA T( AX * ' 6Hl. KOUN T * 5X 5HND1A6*6X 4H1RRG* 6X * AHNPRE * / ) 

8 FORMAT (/55H OUTSIDE D I AME T LR ( i NCHES ) 

1,112.5/) 

9 FORMAT (55H INSIDE DIAMETER! INCHES) 

„ 1 tl2.5/l 

10 FORMAT ( 35H THE ANGLE EXTENDING I HE POCtCt T KLG ION ( UtG. > 

1 E 1 2 • 5 / ) 

11 FORMAT ( 35H IHfc. ANGL t EXTENDING THE LAND Rto I ON ( DEG . ) 

1 E 12 • 5 / ) 

12 FORMAT ( 55H STEP DEPTH ( I NCHES ) 

1 E12.5/J 

13 FORMAT (55+ OUTE- WIDTH OF THE SH-OUD ( I NCHEo > 

+12+5/D 

14 FORMAT ( 55H INNER WIDTH OF THE SH-OUO( INCHES) 

1 E 12 • 5/ ) 

„ 15 FORMAT { 55H CONVERGENCE ERROR 

1 E 1 2 • 5 / ) 

16 FORMAT ( 54H GRID SPAC1NGS IN THE RADIAL D l R lC T ION ( 1 NCHES ) 

17 F QRMAT ( 5 4H GRID SPAC1NGS IN THl C 1 RuUMF lReN T I AL D 1 RECT I ON ( DEG ) 

18 FORMAT ( 55H MINIMUM FILM THICKNESS (INCHES) 

1 E 12 • 5 / ) 

19 FORMAT (55H RE'wOuUT IONS PER SECOND 
1 E 12.5/ ) 

20 0 — IT ( 1H * E;.2.5> 

2 -+2S -E0-S++ 

1 El 2. 5/) 

22 FORMA T ( 55H V I SCOS I T Y ( LB-SEC/ +N4S2 ) 

1 E 12 • 5/ ) 

23 FORMAT ( 55H LA-MQ A ( 6 . *V I $$6 • 28iRPoa> ( RO/HM I N )** 2 /PA ) 

1 E12C5/D 

24 FORMAT ( 1 2H R L3 I I ) * T H ( J ) ) 

25 pO-MAT ( 25H HL » H- * HLH ,HRH * HTH » HBH ) 

2 — 1TU15+ Aiu+T JD TETC 0 

27 F0-MATI45+ F ( ' ) * BU + D * C l +DT A ( +T I - 1 ) , A ( I * l ) * A ( + * ++ 1 ) *10X*2n 

1 ,I5*10X*2HJ=,I5) 

28 FORMAT ( 25H MATRIX IS SINGULAR AT J= l3»i6H,LASt ABANDONED./ i ! - 

29 FORMAT ( 29HOF I NAL PRESSURE U1 sTR i dUT iON. //) 

30 FORMATt //18H CASE CONVERGES TO F9.6, 

1 14H AFTER 13*11H ITER 
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IONS) C11 

31 FORMAT ( / i°( IX * F 1 1 • 7 ) ' 

33 pORMAT£3SH FILM THICKNESS AT 0*D (INCHED) 

1 E12*5/ ) 

34 FORMAT ( /55H TOTAL LOAD(LBS) 

1 E 12 * 5 / ) 

35 FORMAT l 5 5H DIMENSIONLESS L0AD=L0AD/ ( ARE A*PA ) 

1 E 12 « 5 / ) 

36 FORMAT ( 55H HORSEPOWE- LOSS 
1 E 12 • 5 / * 1H 1 ) 

37 FORMAT t 55H STATIC SQUEEZE FILM VEL0 (IN9SEC) 

f 6 T AT +C F + LM V£lOC+TY * + TH HO-.TEP DEPTn 

1 E12.5/) 


NR=5 
NW = 6 

90 READ I NR * 2 ) LAST * NPAD * M #N * I H » JH > I HH • JHH »LKOUNT * NO I AG » I RKG * NPR t 
1 * ND1 

READ t NR * 3 ) DO »D I » THG * THL ♦ STEPD * *0* W I » ERROR 

NN=N~1 

MM=M- 1 

WRITE ( NW * 1 ) 

I F ( IRRG.NE.nGO TO 91 
READ I NR *3 M DRU > t 1 = 1 tMM) 

READ ( NR * 3 > ( DTP { J ) » J= 1 *NN ) 

WRITE(NW*16) 

WRITE! NW. 20 ) ( DR Cl) » 1 =1 .MM ) 

WR I TE ( NW * 1 7 ) 

WRITE (NW* 20 ) < DT H < J ') . J = 1*NN) 

TWD- ( DO-DI ) *0- 5 
IFCNSYM.EQ.H TWD*TWD*0.5 
TlG=THG+THL 
DO 700 1*1 .MM 

700 DR< I )=DR( I )*TWD 
DO 701 J= 1 * NN 

701 DTH(J)-DTH(J)*TLG 
GO TO 95 

91 FMM=M-1 
FNN=N~1 

I HH 1 = I HH- 1 
AIHH1- IHH-IH 


.NSYM 


BIH1*M-IHH 
IHl= IH-1 

aihi=ihi 

IFCNSVM.EQ. 1 ) GO TO 710 
IF ( IH • EQ • 1 ) GO To 709 
DO 92 I = 1 1 1 H i. 

92 DR ( I ) *W I / A I Hi 

709 DO 97 I=1H,IHH1 

9 7 DR( I ) =( CDO-i.*: ' *»5-Wl-WQ ) /AlHHL 
GO TO 720 

710 FIHH1*IHH1 





TEMP* ( (DO-DI )*0*25-W0) /FIHH1 
DO 711 1 3 i ► I HH 1 

711 DR l I ) =TEMP 
720 DO 98 I = I HH *MM 

98 DR(I) s W0/B I HI 

JHH1-JHH-1 
A JHH1 = JHH 1 
BJH=N- JHH 
DO 93 J = l> JHH1 
93 DTH! J ) =THG/AJHH1 
DO 99 J^JHH.NN 

99 DTH(J)»THl/B JH 

95 CONTINUE 

READ ( NR » 4 ) NV I i>M t { V I SA ( l ) * I =1 * NV I SM ) 
READ ( NR * 4 ) NRPE-Mi* ! RPSA! I) .I-l.NRPSM) 
READ ! NR » 4 ) NPAP » { P AA ( I ) . I=1.NPAM> 
READ(NR»4)NHMK. I HMINAt I ) » l*l.NHMM) 
WRITE INPUT DATA 

READ '(NR. 3) CHOUT'C I ) .I-l.NHMM) 

READ ( NR * 4 > NHDTM . ( HDOT ( I) * I - 1 .NHDTM ) 
WRITE(NW,5) 

WRITE!NW,6)LA$T.NPAD*M,Nt IH. JH* IHH.JHH 
WR I TE ( NW * 7 ) 

WR I T t ! NW . 6 ) LK*OUN T .NO I AG. IRRG.NPRt 

WRITECNW.SIDO 

WRITE(NW»9)DI 

WR I TE ( NW» 10 ) TNG 

WRITE (NW.ll ) TH l 

WRITE(NW»12)STEPD 

WRITE(NW.13 )WO 

WR ITE ( NW. 14 ) W I 

WRITE(NW*15) error 

IF(NDIAG) 731 .731.730 

730 WRITE! NW . 16 ) 

WR I TE ( NW . 20 > (DR ( U * 1 = 1. MM) 

WRITE! NW .17) 

WRITE (NW. 20) (DTHTJ) .0*1 .NN) 

731 CONTINUE 

DO 1000 NVIS=1.NVISM 
vl S-VISA! NVIS ) 

DO 1000 NRP$ = 1 . NRPSM 
RPS=RPSA!NRPS) 

DO 1000 NP A* 1 * NPAM 
PA=PAA ( NPA ) 

DO 96 I = 1 *M 
DO 96 J = 1 * N 

96 Q( I »J) =1. 

DO 1000 N+DX31TN+DTM 
DO 1000 NHM-l.NHMM 
HMIN=HMINA(NHM) 

CONE*HOUT!NHNU /HMIN 
CONE 1* CONE- 1 * 

WRITE (NW.33) HOUT(NHM) 
WRITE(NW.18)HM1N 



WRITE(NW,19)RPS 

WRITE(NW,21)PA 

WR I TE ( NW, 22 ) VI S 

0- + TE UNO T 3 L> +DQTUN+DTD 

PI =3 • 1416 

RO= DO/ 2*0 

R I =D I / 2 • 

PLAM=6,*V1 S*RPS*2**P I* ( ( RO/HMIN)**2* I /PA 
CHDOT = HDOT ( NHDT ) / ( 2 . *P I *RPS*HM I N } 
CHDT=HDOT (NHDT ) /( 2 . *P I *RPS*S T EP D I 
WR I TE < NW * 39 ) CHDT 
STEP=STEPD/HM1 N 
WRITE{NW,23 jPLAM 
STEPH=STEP*0. 5 
CC~0. 01745329 
C GENERATE COORDINATES 
T H ( 1 ) - 0 • 

NN=N-1 

DO 100 J - 1 * NN 

100 TH( J+l )=TH( J )+DTH( J)*CC 
R ( 1 ) =DI /DO 
MM=M~ 1 

DO 105 1 = 1, MM 
105 R ( I+l)aR{ I ) +DR ( l >/RO 
D I DO-D 1 /DO 
OD I DO= 1 • — D I DO 
L6=l 
K8 - 1 
ADMY- 1 • 

109 DO 107 J= 1 , N 

107 H(L6*J)=ADMY 

IF ( K8 • EQ • 2 > GO TO 108 
L6 = M 
K8”2 

admy^cone 

GO TO 109 

108 DO 106 1 = 2* MM 
RS9=R{ I ) 

DO 106 J=1,N 

10 6 H( I * J ) =1 ,+CONEl* < ( RS9-DIDO ) /ODI DO ) **2 
iFtNDIAG) 114,114,110 

110 WRITE(NW,24) 

WRITE (NW , 20 ) {R( I 

WRITE (NW, 20 ) (TH(J) ,J=1*M) 

C GENERATE A1(I,J) TO A6M»J) 

DO 112 I ::: 1 * M 

112 WRITE ( NW * 20 ) ( H ( I , J ) * J* 1 ,N > 

114 DO 140 J = 2 * NN 

DZ1=TH( J+D'-THl J-l) 

DZ2=2.*DZ1 

DZ3 = DZ 2* ( TH ( J > -TH ( J-l ) ) 

DZ3= 1 • / DZ 3 

D24~ 1 * / ( DZ 2* < TH< J+l )-THC J) )) 

DO 140 1=1, MM 
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HL=H< I * J-l ) 

HR = H< I * J+l ) 

HLH=(H( 1 .J-l )+HlI ,J) 1 *0 • 5 
HRH=(H( I • J+-1 ) + H ( I * J ) )*Q mb 
HTH= ( H { 1 + 1 * J )+H( I * J ) 1*0*5 
HTEMP=H( 1-1 ♦ J) • 

I F ( I * EQ • 1 ) HTEMP a H { I , J) 

HBH = ( HT EMP + H ( l * J ) ) *0 * 5 

+F ( J- JHH ) 1X1*120*130 
POINTS AT THE LEFT SIOE OF THE STEP 
111 IK(I-IH) 130*1 15*116 
c points on the bottom edge 

115 IFINSYM •EQ. 1 ) GO TO 116 
735 HL-HL+STEPH 

hr-hr+steph 
hlh=hlh+steph 
hRH=HRH+STE PH 
IF ( I.EQ*IH)HTH=HTH+$TEP 
I F ( I.EQ. IHH)HBH=HBH+5TEP 
GO TO 130 

116 IF(I-IHH) 117.735*130 

c points in the pocket 

117 HL-HL+STEP 

hr=hr+step 

HLH=HLH+STEP 

hrh=hrh+step 

hth*hth+step 

HBH=HBH+STEP 
GO TO 130 

120 I F ( +- IH ) 130,125,126 
C BOTTOM OR TOP CORNER 

125 HL=HL+STEPH 

hlh*hlh+$teph 

IF(I.EQ.IH) HTH=HTH+STEPH 
IF(I.EQ.IHH) HBH=HBH+STEPH 
GO TO 130 

126 IF(I-IHH) 127*125*130 

-C POINTS ALONG THE VERTICAL EDGE 

127 HL=HL+STEP 
hlh = hlh+st ep 
hth=hth+steph 

HBH=HBH+$TEPH 
go TO 130 

130 IF(NDIAG) 132,132*131 

131 WR I TE ( NW * 25 ) 

WRITE! NW,20 )HL ,HR*HLH,HRHtHTH*HBH 
C GENERATE A1 ( I * "J ) .ETC 

132 RTEM=R( 1-1) 

DRT*DR ( I“l) 

IF! I .EQ* 1 ) DRT *DR< I ) 

IF t I *EQ.l )RTEM«R( I) 

TEMPl = Pi_AM*R [ l ) /D22 

TEMP2= 1.0 / tA.OMDRl I )>DRT 

Al( I * J ) - TEMP 1 *HL 


)/R0) 
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A2< I * J 3 = TEMP 1#HR - 

A3 ( I * J > =DZ3*HLH**3/R ( I ) 

A4< I »J)=DZ4*HRH**3/R( I ) 

A5 ( I * J ) =TEMP2*HBH**3/ ( DRT /R03 *(R(I)+RT£M J 
A6£ I » J) «TEMP2*HTH**3/( DR ( 1 )/RO) * l R { I) +R ( 1 + 1 J) 

A6 £ I * J ) *PlAM*P ( I J *CHDOT 

1 F { NO I A 6 ) 1 40 ♦ 1 40 . 1 3 3 

133 WRITE(NW,26) 

WRITE(NW*2Q)A1( I * J 3 ♦ A 2 { I > J 3 * A3 ( I > J ) * A4 £ I * J ) * A 5 (I *J)*A6tI*J) « 

1 A8 C I , J ) 

140 CONTINUE 


KOUN T = 1 

141 DO 143 I ~ 1 * M 
G £ I t 1 ) -0 # 

DO 142 K- 1 * M 


142 E t I tKf 1)*0. 

143 CONTINUE 

DO 370 J- 1 * N 
DO 310 I = 1 * M 

lF(J.EQ.1.0R* J.EQ.N) GO Tu 210 
IF( I .EQ. 1 ) GO TO 725 
1 1" £ I ♦Eg. M) GO TO 210 
GO TO 726 

723 I F ( NSYM .NE.13 GO TO 210 
726 SQP=SORTt ABS(0(I »J+1)) ) 

SQM 2 SQR T ( ABS < 0 { I * J- 1 3 > ) 

SOQ-SQR T ( ABS ( Q ( I,JJ ) ) 

B( I )=A3( I » Jl+Al ( I »J)/(2.0*SQM) 
C(I)=A4£I*J) -A 2 ( I * J ) / (2 ,0*3UP ) 

AGUM = A 3 £ I • J 3 +A4 £ I • J J +A 5 { 1 * J 3 + A6 £ I , J j 
F ( I ) -A3 ( I » J ) *0 { I t J- 1 3 +A 1 £ I t J ) *SQM 
OT T 2 Q £ 1-1 ,J) 


IFtI .EQ. 1 .AND. NSYM .EQ. 13 UTT~Q(I+1*J) 

F(l)-F(I)+A 5 l!»J) *QT T -ASUM*Ql I • J ) +A 6 ( 1 * J ) 1 + 1 * J J +A 4 £ It J ) *u 

1 ( I * J + 1 ) -A 2 ( I * J )-*SQP-A 8 ( I »J)*SQQ 
F ( I ) =-F £ I ) 


DO 150 Kal,M 
A ( I * K ) - 0 * 0 

IF (K.EQ. I ) A ( J! * K ) =~ASUM-A 8*0 • 3 /SOQ 
A 6 T T 2 A 6 ( I ♦ J 3 

+ F{+ • EO. 1 .AND. NSYM . EQ. 13 A6 T T-A6 £ I ♦ J ) + A3 < I 
+ 1- (K.EQ.++1 ) A( I ,K ) =A6T T 


IF (K.EQ. I -1 ) A( 1 ,K) =A5 £ I , J ) 


150 CONTINUE 

IF(J.NE.JHH) go TO 305 


J) 


I F ( I *EQ. I H. OR- I • EQ. IHH ) GO TO 151 
IF ( I .GT . IH. AND. I . LT. IHH 3 GO TO 152 
GO TO 305 

151 DMSP-STEP* 0.5 
GO TO 155 

152 DMSP- STEP 

155 SQQ=$QR T ( ABS ( Q ( I * J ) ) ) 

TEMP-A1 ( I . J 3 *DMSP/ ( H ( I * J~1 ) + UMSP 3 
F ( I 3 = F ( I ) -TSMP*SQQ 
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A(I*I)=A(1»I )+TEMP*0«5/SQQ 
30 5 If- (NO I AG) 310*310,306 
306 WRITE! NW *27) I ,J 

WR1TE(NW*2G)K( I ) * U ( I )*C(I ) »All *1 — 1) #A< I *1) *A( 1»I+1) 
GO TO 310 

210 0 ( I ) =0. 

c ( n = 0 . 

FU) =0. 

DO 211 JC«1,M 
A ( I *K ) =0 « 

IF ( I.EO.K)A{ I 

211 CONTINUE 
GO TO 305 

310 CONTINUE 

DO 320 1=1* M 
DO 320 K=1*M 

3 20 QSMA l I *K> = AU *K)+B< I )*E ( 1 *K* J-) 
call mat INV < (JSMA »M*BB*0 *DET » ID ) 

60 TO ( 340*330 ) , ID 
340 DO 360 I = 1 * M 
G< I *J+1 ) =0. 

DO 360 K= i * M 

G( I *J+1 )=G( I *J«fl ) +OSMA ( I * K } * { F (6J-BUC) *G ( K * J ) ) 

E l I J+l ) *~QSMA< I *K)*C( K) 

360 CONTINUE 
370 CONTINUE 
DMA=0. 

DO 380 I = 1 * M 

DMA=AMAX1 (DMA,ABS(G( I *N+1 ) ) ) 

380 DO t I * N ) =G ( I *N+1 ) 

DO 400 J J = 2 > N 
J-N+2- JJ 
DO 400 1=1* M 

dum=o. 

DO 3 0 K=1*M 

390 DUM=DUM+E ( I *K> J )*DO( K* J ) 

DUM=DUM+G( I * J) 

DMA^AMAXl < DMA * ABS ( DUM ) ) 

400 DO ( I * J-l } = D U M 
DO 401 I - 1 * M 
DO 401 J= 1 * N 

401 Q( I ,J)=0( I * J ) +DQ ( I *J) 

+F ( ND+AG ) 40 5 *405* 402 

402 DO 403 J=1*N 

40 3 WRITE(NW,20MD0( I * J) *1 = 1 ,M) 

DO 404 J = 1 * N 

404 wR I TE l NW* 20 ) ( 0 ( I * J I * I = 1 »M ) 

405 CONTINUE 

GO TO 560 
330 WRITE (NW*28 )J 
GO TO 1000 
560 (COUNT = K.OUNT + 1 

I F ( KOUNT •GE.LKOUNT J GO TO 561 
I F ( DMA. GT. ERROR ) GO TO 141 



561 wRlTE<NW*30)ERRORfKOUNT ' 
lF(NPRE*EG.l) WRITE! MW *29) 

DO 576 J = 1 » N 
DO 5 76 1 = 1 9 M 
P( I > J)=SQRT ( ABS ( 0 ( I * J) J ) 

576 CONTINUE 

DO 585 I = 1 * M 

IF ( NPRE.EQ. 1 ) WRI TEINW* 31 >(P(1*J)*J=1.N) 
DO 5 BO J=1#N 

580 QoQt J)=P( I 
NN=N-1 

PP t I) *0* 

DO 581 0=1 * NN 

581 PP ( I) =PP ( I ) +DTH ( J ) * ( OOQC J ) +QuQ( J+l ) ) *0*5 
PP( I ) =R ( I )*CC*PP ( I ) 

IF (ND1 )505»585*8OO 
800 WRITE (NW*3i)PP( I ) 

585 CONTINUE 
FPAD=NP AD 

WLOAD=0. 

DO 589 1=1 » MM 

58 9 WLOAD=WLOAD + DR (I ) *{ PP l I ) +PP ( 1 + 1 ) > *0.5 

WLOAD=WuOAD*RG*PA*FPAD 
IF (NSYM.EO.l )WuOAD=2.0*WLOAD 
WBAR=WL0AD/ t 3 . 1 A 1 6* ( RO* * 2“R I **2 ) *P A J 
IF (N5YM.EQ. 1) WBAR=2.*WBAR 
DO 600 I = 1 » M 

IFU.EQ.IH.OR.I.GT. IH) GO TO 593 
591 FRC ( I) =0. 

nn=n-i 

DO 592 J = 1 » NN 

59 2 FRC( I ) =FRC( I >+DTH( J ) 

FRC< I > =FRC( I ) *• C C * R l 11**3 
GO TO 600 

593 IF tI.GT.IHH) GO TO 591 
FRC ( I ) =0. 

NN=N~1 

DO 594 J= 1 * NN 

IF ( J.EQ. JHH.OR.J.GT. JHH) TEMP=1.0/Hl I , J) 
IF! J.LE.JHH)TEMP«1.0/tH( I » J M-STEP ) 

594 FRCt I)=FRC( n+DTH( J)*TEMP 
FRC I I )=FRC( I )**CC*R( I)#*3 

600 continue 

MM=M-1 

HP0W=0. 

DO 605 I = 1 » MM 
605 HPOW=HPOW+FRC( I J *DR( I ) 

HPOW=HPOW*VIS*6.2832*RP$*-Q**3 
HPOW= j HPQW*FPAD*RPS* 6 0./63GOQ.O 
HPOW-HPOW/HM I N 
I F ( NSYM.EG. 1 ) HPOW=2 .0*HPOW 
WRITEtNW* 34 ) W i.O AD 
WRITE! NWt 35 JW13AR 
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W R I T E ( N W » 3 6 ) HP OW 

1000 CONTINUE 

I F ( LAST ) 90,90,1001 

1001 STOP 
END 


SUB-QUT+NE MATlNVUATNlf BTM1 *L>ETE- *10) . ‘ 

matrix inversion with accompanying SOLUTION OF LINEAR EQUATIONS MAI 10002 

NOVEMBER 1962 S GOOD DAVID TAYLOR MODEL bAS-IN AM MAT 1 MAT^OOOA 


DIMENSION A ( 17 ,17) <*B!17*1)*1NDEX!17*3> 

general form of dimension statement 

EQUIVALENCE ( I ROW * JROW ) , (ICOLU * JCOLU )* (AMAXo I* SWAP) 

INITIALIZATION 

M-Ml 

N-Nl 

10 DETER =1.0 
15 DO 20 J = 1 > N 
20 INDEX (J ,3) = 0 
30 DO 550 1=1, N 

SEARCH FOR PI'-'OT ELEMENT 

AO AMAX=0.0 
45 DO 105 J = 1 * N 

IF ( INDEXtJ.3 )-l ) 60, 105* 60 
60 DO 100 £=1»N 

+ F ( INDEX ( X, 3 )-l ) 60, 100* 715 
60 IF < AM AX -ABS (A<J»K))J 85* 100* 10 0 

5 +-OW3J 
0 +CQLU =n 

AMAX-ABs UA(J,K)) 

100 CONTINUE 
105 CONTINUE 

INDEXt+COLU >3) = INDEXl+COLU *3) +1 

260 INDEX! I *1 ) =IR0W 
270 INDEX! I *2 )= ICOLU 

interchange rows to put pivot element jn diagonal 

130 IF (IROW-ICQlJ ) 140* 310* 140 

140 DETER =-DETER 
150 DO 200 l=1?N 
160 SWAP 5 A ( I ROW » L. ) 

170 A! I ROW * u ) = A ( ICOLU *L) 

200 A! ICOLU » L ) = SWAP 

+ F ( M ) 310* 3 1 0 » 210 
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MAT 10 006 
MA T iOOO 7 
MA l 10009 

MA I 10012 
M A I l 00 1 3 
MA 110014 
MAT 10015 
MAT l 00 i 6 

MAI i 00 1 8 
1*1 A i x 0 0 19 
MAI 10020 
MAT iOGaI 
MA I 10022 
MAT 10023 
m A T 100^4 
MAI 1UU^5 

-I AT 10026 

-i AT j 002 7 

•‘t« ' 1 ’ J 0 2 d 

MA r 100 30 


MAT i 00 3 3 
MAT 10034 

MAI 10036 

VA i 1 003 8 
■~\p. I 1UQ39 
- ' 


A I 1 UOm-3 


. i 0 u. a 
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210 

DO 250 L = 'l* M 

MAT 10048 

2 20 

SWAP-B(+-OW*L) 

MAT 10049 

230 

B t I ROW » l ) =B ( ICOLU 9 L ) 


250 

B ( I COLU 9 L 1 =S WAP 




MAT 10052 


D1V+OE PIVOT ROW BY PIVOT ELtMENT 

MAT 1QQ00 
MAT 1005 3 

310 

PIVOT =AlICOuU » I COLU ) 

DETER=DETER*PIV0T 


330 

A ( + COLU » +COLU } * 1 • 0 


340 

DO 350 u=l»N 

MAT 10057 

350 

A ( I COLU * L ) -A < I COLU »L)/PIVOT 


3 55 

+F(M) 380, 380, 360 

MAT 10059 

360 

DO 370 <_=1*M 

MAT 1006 0 

370 

BUCOLU ,L)=B(IC0LU ,L)/PIVOT 

MAT 10062 


-EDUCE NON-PIVOT -OWS 

MAT 10063 
MAT 10064 

380 

DO 550 l 1 = 1 » N 

MAT 10065 

390 

IFTLI-ICOlU ) 400 » 5 50 , 400 


400 

T = A ( LI * I COLU ) 


420 

A ( L 1 * I COLU } = 0 • 0 
IF(T)430, 550*^30 • 


430 

DO 450 t_ = 1 * N 

MAT 10069 

450 

A(L1»L) =A(L1 ,L )-A{ ICOLU *L)*T 


455 

I F ( M I 550, 5 5 0 » 460 

MAT 10071 

460 

DO 500 i_ = 1 * M 

MAT 10072 

500 

B ( LI »L ) =B ( LI ,L )-b( ICOLU *L>*T 


550 

CONTINUE 

MAT 10074 
MA T 1007 5 


INTERCHANGE CUlUMNS 

MAT 1007b 
MAT 10077 

600 

DO 710 I = 1 * N 

MAT 10078 

610 

L=N+1-I 

MAf 10079 

620 

IF ( INDEX ( L * 1 - I NDEX ( L , 2 ) ) 630* 710* 630 

MAT 10080 

630 

JR0W = INDEX ( L ,1 ) 

MAT 10081 

640 

JCOLU = INDEX (i. ,2) 


650 

DO 705 X=1*N 

MAT 10083 

660 

S W A P = A ( K » JROW ) 

MA T 10064 

670 

A(K» JROW) =A ( K, JCOLU ) 


700 

A (K* JCO»-U )=SV/A^ 


705 

CONTINUE 

MAT 10087 

710 

CONTINUE 

MAT 10088 


DO 7 30 K = 1 ,N 

MAT 10089 


IF ( INDEX! K*3 ) -1) 715,720*715 

MAT 10090 

720 

CONTINUE 

MAT 10093 

730 

CONTINUE 
I D= 1 

■"AT 10094 

740 

-ETU-N 

MAT 10096 

715 

ID =2 

MAT 10091 


GO TO 740 
END 

MAT 10092 



Card 1 Format (72 H) 

Identification Card 
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Card 2 

Last 


NPAD 

M 

N 

Ili 

IHH 

JH 


Format (1615) 

Integer to determine whether additional input data 
are to be read 

Last * 1 p no more input data 

Last » 0 , more input data to be read from state- 

ment 90. 

Number of pads (see Fig.A2) 

Number of grids in the radial direction 

Number of grids in the circumferential direction 

Grid number for the bottom edge of the step (see Fig. A2) 
For NSYM*l f set IH-1 

Grid number for the top edge of the step (see Fig. A2) 

Set JH-1 


JHH 

LKOUNT 

NDIAG 


JJRRG 


NPRE 


Grid number for the left edge of the step 

Maximum number of iterations allowed for the pressure 
calculation (recommended value: 10-20) 

Control for diagnostics 

NDIAG 1 t diagnostics print out 

NDIAG * 0 , no diagnostics 

Control for irregular grids 

IRRG * 1 , read in irregular grid spacings 

IRRG =0 , uniform grid spacing 

Control for printing out pressure profile 

NPRE =1 , print out pressure 

NPRE » 0 , no pressure print out 

This integer is used to control whether pressure 
calculation is made for a full pad or half a pad as 
in the case where the pressure is symetrical about 


NYSM 
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the center line (see Fig.A2). 

NYSM =0 , For calculation covering full pad 

NYSM * 1 , For symmetrical pressure profiles Where 

where calculation is only made for half a pad. 

ND1 - Set ND1 « 0 , for normal runs. 

Card 3 Format (8F10.7) 

DO - Outside diameter of thrust brg., 

d in Fig.Al (in.) 
o 

DX - in Fig.Al (in.) 

THG - Angle extending the pocket region, 8 in Fig.A2, (degrees) 

S 

THL - Angle extending the land region, ^ in Flg.A2, (degrees) 

STEPD -- Depth of the step (in.) 

WO - Outerwidtb of the shreud, in Fig.A2, (in.) 

WI - Innerwidth of the shroud, in Fig.A2 (in.) 

ERROR - Convergence factor for pressure iteration. 

(recommended value: .001 - .0002) 


Card 4 


Format (8F 10.7) 


card is requirt 


DR(I), I ■ 1, MM - Dimensionless irregular gird spaclngs in the radial 

direction. (MM - M~l) 
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Card 5 Format (8F 10.7) 

This card la required only IRRG. 


DIH(J). J « 1, NN - Dimensionless grid spaclngs in the circumferential 

direction. <HN - N-l) 



Fig.A-2 




(O0) 


®L + »G 






+ e. 


Card 6 Format <15, 5X, 7E 10.3) 

NVISH - Total number of viscosities to be investigated in the 
production run. 

VISA(I), 1=1, NVISM 

2 

The arrary of viscosities, < lb- sec/in ) 

Card 7 Format (15, 5X, 7E 10.3) 

URPSM - Total number of angular speeds to be investigated in 

the production run. 

EPSA(I) , I - 1, NRPSM 

The array of angular speeds, (Rev. per sec.) 

Card 8 Format (15, 5X, 7E 10.3) 

HP AM - Total number of ambient pressure to be investigated 

in the production run. 

PAA(I) , 1=1, NPAM 

The arrary of ambient pressure, (PSI) 



Card 9 


Fermat (15, 5X, 7E 10.3) 
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This card reads in the film thickness to be investigated in the 
production run. In the case of a parallel film, the film thickness 
in the land region will be read in the array designated as HMINA(I). 

In the case of non-parallel film with an axisymmetric coning of dlshlrg, 
HMLNA. (I) represents the film thickness at the Inside diameter in the 
land region. The film thickness at the outside diameter is HOUT(I) 
*faicb read in Card 10. 




Fig.A3 


The variables to be read in this card are: 


NHMM 


Total number of film thickness to be investigated in 
the production run. 


HMINA(I) , 1*1, NAMM 


HMENA is the film thickness at the inside diameter in 
the land region (inch) 


Card 10 Format (8F 10.7) 

H0UT(I), 1 -■= 1, NHMM 


HOUT is the film thickness at the outside diameter in 
the land region, (inch) 

Card 11 Format (15, 5X, 7E 10.4/(8E10.4) 

N* 113 ™ * Total number of velocity or time variation of gas 

film thickness to be investigated in the production run. 

HDOT(I) - Array of velocity or time variation of gas film thickness 
in in. /sec. 



ro rv> f\; jv 
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PROGRAM RSGALN( INPUT .OUTPUT) 

DIMENSION WN(30I* QN<7), AG(9)*A0G<9) 

D+MENS+ON AHAdOOD. FlA(lOO). F2AU00). F3A(100)* F4AU0G). F 1 1 9 ) 
DIMENSION FbA(100)f F6 A ( 100 ) • F 7 A ( 100 ) . F8AU00I* F9AC100) 

THE TABLES ARE PREPARED IN THE FOLLOWING ORDER THAT l F J A ( I > * J** 1 * 9 ) 

= ( D2*I1A* Dl*14* 02*I4A* -U2*I1A0* -U2*I4AU* U2* 1 5 /HO**2 • 5 

. DI*I5A/HO**2.5 . -DZ*I5AO/HO**2.5 >• 

DATA F1A( 1 ) /O./ *F2A( 1 ) /7.854Z *F3A<1 ) / 3 . 14 16/ *F4A U ) /O. / 

DATA ( F 1A ( I) .1=2.60/ 7 * 8 5 5 5 28 1 E-0 2 . 1* 5 72034 1E-Q1 * 2 * 3603 7 50E-Q 1 • 


2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 

2 


3.15151236-01 
6.36311466-01 
9.6978133E-01 
1*32247156+00 
1.7022791E+00 
2*11856096+00 
2 6 1 1 E+OOT 
3*10 27486+00 
3*71667666+00 
4*42 52U6E+C0 
5*28105296+00 
6*317 42036+0 0 
7. 60441756+00 
9*23868746+00 


DATAUF 1 AU+ ) .1=61.100)/ 


1* 1367114 E + 0 1 
1 *42220386+01 
1.81890236+01 
2 • 394670 7E+0 1 
3.2781965E+01 
4.73852986+01 
7*41901 5'76+C) 1 
1.31 5000E+02 
2. 96569976+02 
1.1832960E+03T 


3 * 94639 126-01 
7*18268016-01 
1 * 05 58 636E + 00 
1*41454486+00 

1 * 8025243E + 00 

2 * 2296 1 466 + 00 
2C707 04 6+OOT 
3. 2527020E+OU 
3C 8387 16+00 
4C6278602E+00 
5*520648 3 E +00 
6.61253726+00 
7*975 746 5E + 00 
9*71721836+00 


1.2001018E+01 
1 • 5 08957 3E + U 1 
1*942430 6+01 
2.57954656+01 
3*5 73462 7E+01 
5*25 3947 36 + 01 
8.44332126+01 
1.57010 56+02 
3*871653 E+02 
2.10183026+03 


DATA( F2A(I ) .1=2.60) / 7.058621 


7*928 49 7 5E+ 00 
8. 1556360E+00 
8*54649766+00 
9. 12072306+00 
9*90838436+00 
1*09530376+01 
1*23164536+01 
1C4085 56+01 
1*63854076+01 
1.9393624E+31 
2.33725636+01 
2*87154976+01 
3*602 3656+01 
4.6282353E+D1 


DAT A ( F2 A ( I ) . I =6 1 T 1 00 ) / 


6*10808566+01 

8.32107376+01 

1*178^1276+02 


7 ■ 97067 3 7E + 00 
8*237 39 22E + 00 
8.67188596+00 
9.29618716+00 
1.01433 86+01 
1.12611046+01 
1 * 27 166406 + U 1 
1C4605207E+01 
1.70624936+01 
2.02847126+01 
2.45612006+01 
3.03290686+01 
3.82681 50E+0 1 
4 * 947 2941 E+U 1 


6.57794356+01 

9.04233876+01 

1.29461606+02 


4.74596806-01 
8.01094116-01 
1.14325846+00 
I. 50844 156+00 
1*905203 1E+0U 
2.34385296+00 
2 CB 3 7 1 83 56 + 00 
3.40149646+00 
4C0580843E+00 
4. 3543216+00 
5.772583 7E + 00 
6 • 924447 66+00 
8. 37044916+00 
1.02291876+01 

1.26843966+01 
1 * 60 334 1 2E+0 1 
2.07835806+01 
2.78596296+01 
3.90962876+01 
5.05720956+01 
9.69408266+01 
i. 89941016+02 
5.26667446+02 
4.7245976E+03 
E + QO * 7.8 72555 
0.02247956+00 
0 . 3295914E+0C 
8 .80908396+00 
9.48551356+0 0 
1.03952236+01 
1.15901526+01 
i.3l43b996+01 
1C5 1597766+01 
1 . 7786946E+01 
2.12409956+01 
2.58419106+01 
3.20764056+01 
4.070761 3E+01 
5.29759806+01 

7.09863796+01 
.85004186+01 
1 .42b0 12 36 + 02 


5.55121336-01* 

0.8 48V 5 6 06- 0 1 * 
1.23208616+00. 

1 • 60430 34E+00 * 
2.01048796+00* 
2.46149196+00* 

2. 7087456+00* 
3.55602376+00* 

4.23 78766+00* 
5.0529163E+00* 
6.03782246+00* 

7. 2545570 £+00* 
8.79064296+00* 
1.07779086+01/ 

1.34226736+01. 
1.70629396+01 * 
2.2284206E+01* 
3.01742856+01* 

4. 29465306+01* 

6. 56954626 + 01 * 

1 . 1 243492^6+0 2 * 
2.34421576+02* 
7.57889766+02* 
1.88783816+04/ 
16+00. 7 . 8958 2 39 E +00* 
8.08407126+00* 
8.43251936+00* 
8.95853386+00* 
9.68934516+00* 
1.06647726+01* 
1.19414676+01* 
1.35994536+01* 
1.573221 £+01* 
1.85625866+01 » 
2.22681706+01* 
2.72234226+01* 
3.39713156+01* 

4. 3370409E+01* 
5.68303396+01/ 

7*67720246+01* 

1 * 07505686 + 0 ^ » 

1 . 5 f / * 



ORIGINAL PAGE IS 
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1 1.7bl<!021E + 02. 1 .9511517E + U2 . 2 . 1B299V 9E+02 . ' 2 .4&33690E+02 * 

2 2 • 770607 16+02 * 3.14332 66+02* 3 • 39 1 1 7 3 16+Q2 * 4 « 1 2585996+02 * 

2 4 • 77 267 36E+02 * 5.56263826+02* 6 . 6 3767 7 36+02 * 7.75536246+02* 

2 9.29617356+02* 1.12748 56+03. 1.38589386+03* 1 • 7 2988 2 76+0 3 * 

2 Z* 19806 106+03* 2.85197376+03* 3. 793728 36+03 * 5 • 20095646+0 3 * 

2 7C400 3266+03% 1C 1 04 1 00 1 6 + 04T 1C 75224476 + 04 * 3 • 0 26 10246+04 * 

2 5 • 906 86 84 E +04 * 1 . 3 99 3 1 3 5 £ + 05 » 4 . 7 19842 36 + 0 5 * 3.77357286+06/ 

+ + OF A0 + + 0+ 2Q60D/ 3 . 1 4 19 5 i Ofc + 00 * 3 . 1 42 98246 + 00 * 3 • 1 4470 2 7E+Q0 , 
2 3. 1471 142E + 00, 3 . 1 50 2 1 98E + 00 * 3 . 1 5402 3 56 + 00 » 3. 15852986+00 * 

Z 3.163/4466+00, 3.16967436+00* 3 . 1 76 3 26 36+QG * 3 • 1 83 7Q9 66 + 00 * 

2 3.19183316+00* 3 . 2 00 70 7 5£+00 * 3 . 2 1 034426 + 00 * 3 • 22075 58E+00 * 

2 3.23193626 + 00* 3 • 24396 02 1 + 00 * 3 . 2367839E+00 * 3*2 7044486+00* 

2 3. 284V616E+00 * 3.30035466 + 00* 3 . 3 1 6'643 36+00 * 3 • 33385696+00 * 

2 3.35201446+00* 3.3/114426+00* 3.3912749E+00* 3.41243696+00* 

2 3+4346625E+C0* 3.45/98656+00* 3 . 4b 2445 86+00 * 3 • 50808U 1 E+QQ * 

2 3.53493136+00* 3.36304516+00* 3.59246916+00* 3 • 6 2 32 549 E+00 * 

2 3.65545766+00* 3.68913606+00* 3 . 7 2435 3 1E+00 * 3 . 76 1 1 76 3 6+00 t 

2 3*79967786+00* 3.8 3993 326+00* 3 . 88 203 1 3 6 + 00 * 3 . 926036 1E + 00 * 

2 3*9721 049E+00* 4 . 020 2 8 1 2 t + 0Q » 4.070696 36 + 00* 4. 1 234 70 2E + 00 * 

2 4.17873236+00 j 4 . 2 366 2 2 8 t+00 * 4 , 29 729 3 26+00 > 4. 36 090 78 E+00 * 

2 4.42 764496+0 0* 4 . 497 69 846 + 00 * 4. 3 7 12 794E + 00 ♦ 4. 6486 1 8 2 E + 00 ♦ 

2 4.72996636+00* 4.81559 It+OQ* 4. 9058 18 8E+00 * 5 . 0009572E+00/ 

DAT A ( F 3A ( I ) , I =61 *100)/ 

2 3 . 101 38u 16+00* 5.20749116+00* 3.31973706+00* 3 • 43861 34E+00 * 

2 5.564b716E + C0* 5.6985 2636+00* 3 . 8408664E + C0 > 5.99246346+00* 

2 6.1341 7 56+00* 6 . 32 7 0 2 1 5 1 +00 * 8.31206066+00* 6. 7 1Q6352E + QQ* 

2 6.92413826 + 00* 7.1542 5996 + 00* / • 4029 28 8E+00 > 7.6 7238486 + 00* 

2 7 .96524296+00* 8 . 2845 764E+UQ . 8 . 6 3402046+00 * 9.01790696+00* 

2 .44143 06 + 00* .91092126 + 00* 1 . 0434Q6 7E + 0 1 * 1. 10204 1 3E+0 1 » 

2 1. 1681892E + 01 * 1.243363 36+01 * 1 . 329 3 1 i 4£ + 0 i , 1 . 42918 i 8E+0 1 * 

2 1 .5457768 E+ 01* 1. 6839404 E+ 01* I . 8501 808E+0 1 * 2 . 0 5390 8 9E + 0 1 * 

2 2.30925606 + 01 , 2.63845 286 + 01 * 3 . 0 7859 1 4E + 0 1 * 3 . 696 5 1 3 7E + 0 1 * 

2 4.626058 3 E + 01 * 6 1 7988926 + 0 1 » 9 . 29 7 1 744E + 0 1 * 1 . 868 02 1 3 E + 0 2 / 

+ AT1UF4AU + DT + 32 T60D/ 6. 736313E-02* 1.3753 146-01* 2 . 06 5 50 24E-0 1 * 
2 2.737987 IE- 01 , 3*45390216-01* 4 . 1 54 1 249E-0 1 * 4.859 545 86-01* 

2 3 . 57 1U716E-01 . 6 . 2 896 2 8 06-0 1 » 7.01616236-01* 7 . 75 1648 1 E- 0 1 * 

2 8.49708636-01 * 9 • 2535 1 026-0 1 ♦ 1.00219896 + 00* I • QBO36316+Q0 » 

2 1.15995896+00, 1 • 241 1 062 t+00 ♦ 1.32393046+00, 1.40856266+00* 

2 1.495 14016 + 00* 1.58 380 7 36 + 00* 1 . 6 747 1 6 1E + 00 » 1 . 76802656+00 * 

2 1.86390776+00, 1.96 2 5 3856 + 00* 2 . 064108 56+00 * 2. 1688 189E+00 * 

2 2.27688356 + 00, 2.3885 3006+00* 2 . 50400 1 4E+0U ♦ 2. 6235 5706+00 * 

2 2.7474742E+00, 2 • 8 760 5 026 + 00 * 3 .0096038E + 00 * 3.14847 766 + 00* 

2 3.29304006 + 00* 3. 44368 8 2 E + 00 , 3 . 600850 7E+00* 3.764990 76 + 00* 

2 3.936 60956+30* 4. 1 16 25046+00 , 4 . 304503 7E+00 * 4.50201106+00* 

2 4.709471 1E+00* 4. 92 7646 7E+U0 , 5 . 1 57 3 7 1 3E+ 00 * 5.3995 5 7 36+00* 

2 5.65520546+00* 5.92341546 + 00* o . 2 i 1397 86+00 • 6. 5i44o79E + Q0* 

2 6.83616206 + 00* 7. 17805486 + 00* 7 . 34 1 98 1 2 E + 00 * 7 . 9299605E+00 * 

2 8.344244 /t+00* 8. /8 735136 + 00* 9 . 262 102 4t+00 * 9. 7 7166946 + 00/ 

DAT A ( F 4A { I) * 1=61 *100) / 

2 IC031 6276+01 * 1.0 100136 + oiT 1 C 1 5474 i 3E + 0 i * i . 2 23 7038E+Q 1 * 

2 1.2984832E+01 * 1.37976116+01* 1.46832106+01* 1.56306736+01* 

2 1.67 104876+01, 1 • /8 748636+0 1 * 1.913809 16+01* 2 . 03769 796+0 1 * 

2 2.21513986+01* 2 • 39049 866+0 1 * 2.38660296+01* 2.80686066+01* 

2 3C03540 /06+ 01 , 3C33/2 16 + OlT 3C658 7 7 3 1 6 + 0 1 * 4.0273 056 + 01 , 
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2 4.45350706+01, 4*94892386+01 » 5 * 5 2 978 7 2E + 0 1 * 6 • 2 1688 7 3E + 0 1 * 

2 7.03772696+01 , 8*02931 0 7c + 0 A , 9 . 2 42342 2 E + 0 1 * 1.07478196+02, 

2 1.264/5906+02* 1 . 309 1 9476+02 * 1 . 8 3 1 0 1 2 0E + 02 . 2 . 2665762E+Q2 » 

2 2 ♦ 8 7 6 3 8 3 7 6 +0 2 * 3 • 76 7 1 5686+02 * 5 • 1 4 1 69 9 5E+0 2 » 7 • 424904 1 E + 02 * 

2 1 • 163488 16 + 03 * 2 .07457536+03 * 4 . 6822 3 3 3E+03 * 1 » 8789967E+04/ 

DAT A ( F 5 A { I ) *1=1* 39)/ 0* » 2 • 749 /442E-Q 1 » 3 • 3045962 E-Q 1 * 

2 8. 2696846E-01 * 1 . 1 05 0 1 80E+00 * 1 • 385 1 3 1 5E + 00 * 1 • 6678408E+00 * 

2 1.95368836+00* 2 . 2 4 3 2 2 9 3 E + 00 * ^ ♦ 3 3703 3 3E + 00 * 2 * 8 356944E+00 * 

2 3.13 816 1 E+00 * 3 . 45UQ 3 30 6 + 00 * 3C 767004 1E + 00 , 4 • 09 14 1 7 8E+00 * 

2 4*42399336 + 0 0 * 4 . 76549 45 E + 00 * 3 . 1 16 7 1 3 6E + 0U * 5 . 4784962E + 00 * 

2 5.8517341E+00, 6 . 2373745E + 00 , 6 • 6364245E+00 , 7 • Q499568E+Q0 * 

2 7.47911716+00* 7.92513126+00, 8 . 3893 1 2 5E+00 * 8 . 87307 16E+00 , 

2 9.3779257E+00* 9* 9055097E + 00 * 1 .0457588E+01 1 • 1036068E+01 * 

2 1C1643017E+01 * 1 C 2 2 806 7 3 E +0 i T 1C295 147 3E+0 IT 1 . 3658065E+Q 1 * 

2 1 .44Q3336E + 01 , 1 . 5 1 90438 E + 0 1 * i • 6022 8 1 8E+ 0 1 * 1 . 69042 50 E+0 1 » 

2 1C 7 38 776+01* J.C 831251t + 0iT 1C 886382E + 01T 2 • 1009798E+0 1 * 

2 2* 2207604E+01 * 2*34865606 + 01 * 2 . 4854161E + 0 1 * 2*63187386 + 01 * 

2 2*7889564 1 +0 1 * 2*95/69906+01, 3.13925606+01, 3 • 3 349 52 2 E+0 1 * 

2 3. 34617636 + 01 » 3 . / 7463 1 36 + 01 , + • 02 2 148 1>6 + 0 1 * 4* 2 9082 34E + Q i • 

2 4.58303326+01 , 4.90148 366 + 01 , 3 * 249262 2 E+0 1 » 5 • 62990 36E+0 1 / 

DATA < F5 A ( I ) , f»60 * 100 ) / 

2 6.047 46 7 1 E+0 1 * 6.50662986 + 01 , / » 0 1 2 6000E + 0 1 * 7 • 5 722 544E + Q 1 » 

2 8.19230586 + 01 * 8.8815087E + 01, V . 6499 1 3 3E+0 1 » 1.05Q938QE + G2 , 

2 1* 14739766 + 02 * 1 * 2 560462 6 + 02 * i • 3 7889 2 8E + 02 * 1 • 5 18358 1E+02 * 

2 1.677377 5 E+0 2 * 1 * 85953326 + 02 * 2 . 0692 2 74E + 02 * 2 • 3 1 19 1 03E+02 * 

2 2*59438386+02, 2.92521056+02* 3 . 3 1 52694E+0 2 » 3. 7785238E+G2 * 

2 4. 33309266+02 * 5.00277146 + 02, 5 • 8 192 2356 + 02 * 6 • 8 25 19 7 8E+02 * 

2 8.0793374E+02* 9.66354146 + 02* 1 . 169449 1 E + 0 3 * 1 • 434220 7E + 03 , 

2 1,78608526+03* 2 . 264 168 7 6 + 0 3 , 2 . 9307 839E + 0 3 * 3. 8892080E+03 * 

2 5 . 31892426 + 03 * 7 . 55028296 + 03 * 1 • 1236044E + Q4 , 1 . 7 78 7 769E + 04 * 

2 3*06426976 + 04* 5*966417 7 6 + 04* i . 4098 7 9 8E + C 5 * 4. 7435606E + 05 , 

2 3*78303316+06/ 

DATA ( F6A ( I ) * ( = 1 * 59)/ v . 8 540 1 8 <+ E+00 » 7 • 85 5 5648E + 00 * 7 . 86Q2Q 7 1 E+00 * 
2 7.8679533E+00* 7.87881756+00* 7 • 8928 1 9 1E+Q0 » 7 . 90998 34E+00 • 

2 7.9303414E+00* 7.95393006+00* 7 . 980 792 4E+00 * 8 . 0 109 77 9E + 00 , 

2 8.0443422E+00* 8C08 1 5478E + 00 * , 1 2 2 0643E+00 * Q • 1 66 1662E + 0Q » 

2 8. 2139440E + 00* 8.265483 E + 00* 8 • 320888 5E + 00 , 8 • 3802672E+00 * 

2 8.4437389E+00 * 8 • 5 1 1432 1 E+00 * 6 . 58 3486 OE+OO * 8 .6600507E+00 » 

2 8C7412 3E+0C* C 273737E+00T C 184 52E+00* . 0 148 55 5E+00 * 

2 9.1 1667326 + 1)0* 9 . 2241 836E + 0Q * 9 . 3376399E + 0G * 9 • 4573 1 50E+Q0 * 

2 *583502 6+00* C7 1 652 05E+00 * * 567093E + Q0* 1 • 0004438E+0 1 * 

2 1.01601 046 + 01 * 1 .0324138E + 01 , 1C0497005E+0 1 * 1 • 0679206E+0 1 • 

2 1.087 1287E + 01 , 1 . 1073838E + 01 , 1 . 128750 IE+OV* 1 . 1512970^+01 * 

2 1. 1751005E+01 , 1. 20024306 + 01 , i . 2268 144E + 0 1 » 1 • 2 549 1 326+Q 1 , 

2 1.2 4646 7 E + 01 * 1C316132 E + 01* 1C349501 16 + 01 * 1 • 3848932E+01 ♦ 

2 1.42246586+01 * 1 • 46239 16t + 0 1 * i • 50486 1 5E+0 1 » 1 . 5 50086 8E+0 1 * 

2 1 .59830246+01 * 1.64976 3E + 01 * 1 . 7047788E + 01 * 1 . 7636566E+Q 1 / 

DATA ( F6A{ I ) ,1 =60,100)/ 

1 » 8 2 6 76 7 HE +01 , 1 • 8 9 4 5 2 2 6E +0 1 , 1 • 9673837E + 0 1* 2 • 04587406+0 1 » 
2.13058676+01* 2 . 222 1 97 1 E + 0 1 * 2 * 32 i 4 764E + 0 1 * 2 . 429308 5E+0 1 ♦ 

2. * 546 7 1 09 6 + 01 * 2*6 7486006+01, ^ . 8 1 5 1 209E + 0 1 * 2 * 9690866E+Q 1 » 
3.138624 3 6+01* 3 . 32 59 3 5 2 6+0 1 * 3 . 5 3 362 90E + 0 1 * 3. 7648184E + 01 * 

2 4C023241 76+01 * 4C 3 1 342 01 E + 0 1 T 4C640 644E+01T 5 • 0 1 2348 1E+Q 1 * 

2 5.4 36 27 336 + 01 * 5*92316606+01, O * 466 6+0 1 * 7 . 142942 1 6+0 i » 
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2 7.915U96E + G1* 8.832 165 3t+Ul * 9 • 9 3 33 2 28t+Q 1 * ' 1 • T272 193E+Q2 * 

2 1.2 235 5 / E + 02 * 1C4 4 32 1 E + U2 , 1C 76416b 3E + 0 2 t 2. 1 1Q456GE+02 * 

2 2.5760616E+02* 3 • 2 2 3 58 70 E + 02 * 4 . 1630690E+02 * 5 . 6028454E+02 • 

2 7.97778746+02* 1 . 232600 lfc + 03 * 2 . 1666353E + 03 * 4. 821018QE+Q3 • 

2 1. 06 0 72E+04 / 

DATAUF7AU+) * + = 1 * 5 )/ 6C2 32G00E+00* b. 2845747E+GQ*6*2887012E+00* 


2 6.2955873E+00* 6 . 305 2457E+U0 

2 6.35196 1 8E+00 * 6,37204316+00 

2 6.45 2 67 27E+00 * 6 . 4856 166E+00 

2 6.60375146+00, 6Cb4 4 lit+QO 

2 6. 808 50 18 E+00* 6 . 8689366E+00 

2 7 .0743452 E+00 , 7 • 1 5 1 3 5 19E+O0 

2 7.410440 1E+00. 7 * 506Q 376E+00 

2 7. 82 9448 2 E+00* 7 . 94904 0 1 E+00 

2 8. 348462 3 E + 0 0 » 8 • 496 3949 E +00 

2 8. 990 684 2 E + 0 0 . 9 • 1 7 39 8 4 7 1 + 00 

2 9. 78 7 96 2 2 E+00* 1 . 00 1 6 362 E+0 1 

2 1.078 47 18E + 01, 1.1071 4 £ + 01 

2 1. 2044 190 E+0 1, 1 . 241 0079E+0 1 

2 1.3658790E+01* 1.4132695E+01 

DATA ( F7A{ I ) , I ==6 J * 100 ) / 

2 1.5767185E+01 . 1 . 6 39452 1E+0 1 

2 1.8586546E+01* 1 . 94396 2 1 E+0 1 

2 2 . 247 u86 3 E + 0 1 , 2.36/14 lt + OI 

2 2. 039 Llh+01* 2. 7 726 E+Ol 

2 3.6 3962 iOE + Ol * 3. 9131565 E + Ol 

2 4C 53764E+01T 5C4510 25E+01 
2 7 • 364097 5 1 + 0 1 * 6 . 2480692E + 01 

2 1.2208965E+02* 1 . 42 2 1 4 3 1 E+02 

2 2.47 3 3660E + 02 * 3 . 1 08 1 241 E + 02 

2 7C7 2 615E+02T 1C20 46 E+03 

2 1.897567 1E+34/ 

D AT A ( F6 A ( I) * 1 = 1 * 59 )/ 0. 

2 8.26 6846 E — 31 * 1 . 1 05 0 1 80E + 00 

2 1.9536883E+00, 2 . 243229 5 t+00 

2 3.13 8 16 1E+30 * 3, 450 03 30E + U0 

2 4.423 9933E+J0* < 4.7654945E+U0 

2 5.851 /341E+J0* 6 • 2 3 7 3 74 5 E+ UO 

2 7 .47911 7 1E+00 * 7 . 9 2 5 1 3 1 2 t+00 

2 9. 377925 7E+ JO* 9 . 905 5 09 7E + 00 

2 1. 16430 17E+01 * 1. 22806 73E+01 

2 1 .440 3 3 36 E+0 1 , 1 . 5 190438E+01 

2 1 .7838877E+01 * 1 . 8 63 1 2 5 1 E+0 1 

2 2.2 2076O4E+0I, 2C 3486 5 6UE+0 1 

2 2.7889564E+01* 2.95/6990E+01 

2 3.5461763E+01* i.7/463l3t+Ul 

2 4 • 58 30 3 3 2 E+0 1 > 4 . 90 1 4 8 36E+0 1 

DA TA(F8A( + ) *1=60*100)/ 

2 6.047467 1E+01 * b . 5 066298 E+0 1 

2 8. 1V23068E+U1 » 8 . 66 1 5 08 7 t+U 1 

2 1.14 7 3 9 / 6 E + 0 2 * 1.256Q462E+02 

2 1*67737/56+02* 1.6595 3 3 2 E + u 2 

2 2. 5 943 8 38 E 4 0,' , 2.9262105E+O2 

2 4.3330 26t4iv, 3 * 002 77 1 4t + U2 


6 *3176945 E+00* 6. 33295 70E+00 . 

6 . 395940 6 E + 00 * 6 . 4227999E+00 * 
b.5216958E+00, 6 • 56098 1 5E + G0 t 
6.6988937E+00*; 6 • 7 5 16606E+00 » 

6 .933294 0E+ 00* 7 .00 171 34E+00 * 

7 .2 32908 4E+00 * 7 . 3 192034E+Q0 * 

7 . 60863 1 9E + 00 • ' 7 . 7 I60774E+00 • 
8.0751721E+00* 8 . 208 1888E+00 * 
6.6524219E+00* 8 . 8 17 0 1 47E+00 * 
9.367518 4E + 0Q * 9 . 5 7 19399E+Q0 * 

1 .02 5 796 9E + 0 1 * 1 . 0 5 1 3769 E+0 1 * 

1.13 7 668 5 E + 0 1 * 1 . 1 7002 7 1E + 0 1 * 
1.2799 7 51E + 01* 1 . 3 2 152 14E+0 1 * 

1 .48399746+01 . 1.5183642E+01/ 

1 .7 07 00 7bE + 0i * 1 • 7 7986D6E+Q 1 * 
2.0365486E+01* 2 . 137263 7E+Q 1 , 

2.49b 7684E + 01 » 2 . 64348 1 1 E + 0 1 * 
3C1759145E + 01 * 3 • 3946705E+Q1 * 

4 .22b6341£+oi* 4. 58064 5 6 E+0 1 * 

5. 90 lllE+OiT 6.6212373E + 01 • 
9.3116396E+01* 1 .0607436E + 02 . 

1 .6799712E + 02* 2 . 0 179468E+02 • 
4*031146 3 E+0 2 » 5 • 4489 1 5 7E+02 » 
2C1359358E+03T 4.7745 32 1E+03* 

* 2.7497442E-01*5.5045962E-01» 
1. 388131 5E+00* 1. 667840 8E+00* 

2 • 5 3703 5 3E+00 * 2 . 8356944E+00 * 
5.7670041E+00* 4.0914178E+00* 
5.11671 38E + 00 * 5 . 478496 2E+00 * 
o. 636424 5E+0U* 7 . 049956 8E+00 * 
8.3893125E+00* 6 • 8 730 7 16E+00 * 
1.04575686+01* 1. 1036U68E+01 » 
1.2951473E + 0 1* 1 . 365806 5 E+0 1 * 

1 .60228166 + 01 * 1 .6904250E+01 * 
1.9886 3 8 2E + 01* 2 . 1009798E+0 1 * 
2.4854161E + 01* 2 . 63 18 738E+0 1 * 

3. 1392 5 86E+01. 3 . 3349322E+0 1 * 

4.02 214866 + 01* 4 • 29082 3 4E+0 1 ♦ 
5.2492622E+01* 5 . 6299036E+0 1 / 

7.0 1 28000E + 0 1 » 7 * 5 7 2 2 544 E + 0 1 » 
9.6499 133t+0l * 1 • 0509380E+02 » 

1.3 7889 2 8E + 0 2* 1 . 5 18358 1E + 02 * 
2.0692 2 / 46 + 02 • 2 . 3 1 1 9 1 0 3E + 02 * 
3.31526946+02* 3 • 77852 3 8E+02 * 

.8 » 92^.3 66 + 02 » 6.8251978E+U2* 
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2 8.07 3374E + 02 * 9.6635414E+02 * 1 • 1694491E+03* 1 . 434220 7E+03 • 

2 1. 78608b2E+03, 2 . 264 1 68 7 E + 0 3 . 2 . 9 30 78 39E+0 3 * 3 . 8892080E+03 * 

2 b. 3189242E + 03* 7 . b 50 2 8 29E + 0 3 « 1 . 1 236044E+04 , 1 . 7 78 7 769E+04 , 

2 3.0642697t+04, b • 9664 1 7 7 1 + 04 * i • 4098 7 98E+0b » 4. 74356Q6E+0 5 ♦ 

2 3.7830331E+06/ 

OATAUF AU + DT + 3 1 T b U/ 107 0 000E+01# i . b 7 1 4 186E+G i 1 1 . b 732 762 E+0 1 * 
2 1 • 57637 77E + 01 , 1 . b 80 7 3 1 b t + 0 1 « 1 • 5863493E+G 1 * ' 1 • 5932463E+01 * 

2 1 .60 1441 3E + 01 * 1 .6 lU9566t+ui * i .62 1 8 1 8 3E + 0 1 * 1 • 6340569E+Q 1 # 

2 1.6477061E + 01 * 1 . 66 2 8 04b E+ 0 1 • 1 *6 7 939 bOE + O 1 » 1 • 6975 2 5 2 E+0 1 • 

2 1. 71724 78E+01 * 1 • 7 3862 0 7 L + U 1 • 1 • 7bl 7076E + 0 1 » 1 • 7 865 78 1 E + 0 1 * 

2 1 .8133084E+01 * 1 * 84 1 98 1 6E + 0 1 * 1 • 8726884E+0 1 * 1 • 9Q552 74E + 0 1 * 

2 1 .9406U6UE+U1 > 1 * 97804 1 1 E+0 1 * 2 • 0 1 79b99E + Q 1 » 2 . 060b007E+0 1 * 

2 2. 1Q5814QE+01 * 2 * 1 3406 3 7E + 0 1 , 2 . 2054280E + 0 1 • 2 • 260 10 14E + 0 1 * 

2 2. 31829b6E + 01 * 2 . , 3802416E+01 * 2 .446 19l£E+0 1 * 2 . b 1642 14E+0 1 * 

2 2.5912323E + 01 • 2 * 6 709545E+U 1 * 2 . 7559497E+0 1 » 2 . 84661 52E+Q 1 , 

2 2.9433873E+01 * 3 . 046 746 2 t +0 1 » 3 . 1 b 7 22 1 3E+0 1 > 3 • 2 75 3965 t+Q 1 # 

2 3.40191 75E+01 * 3 . 5374992 E + 0 1 * 3 . 6829344E + 0 1 * 3 * 8 39 104 1 £+0 1 * 

2 4.0069889E+01 » <+♦ 18 '/6827E+01 t 4 . 382 40 7 9t + 0 1 * 4. 5925337E+Q1 » 

2 4. 819b974E+0l * b ♦ 06b 32 8 1 E+O 1 * b • 33 1676 bE + O 1 * 5 • 6 20848 1 E+0 1 * 

2 5.9353433E + 01 , 6 . 27800 45E + U 1 * 6.652Q729E+01 * 7.0612547E+Q1 / 

DATA ( F9A ( I ) *1-60*100)/ 

2 7.5098G17E+01 * 8 . 0026082E + 0 1 * 8 . 645 32 7 1E + 0 1 * 9 . 144b 1 1 9E + 0 1 » 

2 9.8077892E + 01 * 1 • 054407 1 E + 02 * i . 1 3638 1 bE+O 2 * 1 • 2 279 3 50E+Q 2 * 

2 1 • 3305279E+02 * 1 .44b8987E + 02 * 1 • b 76 1 28 1 E + 02 * 1* 7237209E+02 ♦ 

2 1C8 171Q8E+02T 2 * 0837956E + 02 * 2 . 3043 146E+02* 2* 35948 12E+02 * 

2 2 • 855693 1E + 32 , 3 • 2 0 1 949 1 E+ 02 * 3 . 6094 1 6QE + 0 2 * 4 . 0 9241 00 E+0 2 • 

2 4.66 48 7 3 E+02 * b . 3649902 E + 02 » 6 • 2 1 1 2 7 40E+02 * 7* 25 197 15E + 02 * 

2 8. 5468744E + 02 * 1 *01793 2E + 03* i • 226822 7E + Q 3 * 1 .4986 157E + 03 * 

2 1.8591183E+03* 2 . 3480042E+03 > 3 . 028 390 5E+ 03 * 4. QQ47739E + Q3 * 

2 b • 4b8 b62 6E + 0 3 » 7 • / 2 3 29 1 3 1 + 0 3 * j. • 145 7 3Q8t + Q4 * 1 • 8Q8 2 7 32 E + Q4 * 

2 3 • 1 0 5 8 8 0 .3 E + 0 4 * 6 . 0 30 1 2 84E + 04 * 1 .42098 1 8E + 03 . 4 . 768032 7E+05 » 

2 3.792641 9E+06 / 

10 FORMAT ( 9 ( E 1 4 .» 7 ) */J 

11 FORMAT (415) 

12 FORMAT l 8 F 1 0 . ) 

13 FO-MAT ( 54H THE AMPL+TUDE OF -ES-ONSE = 

1 E 14*7 * / ) 

14 FORMAT (54H PHASE ANGLE DlFFtRENCE (OtGRbt) = 

1 E14.7 */// ) 

15 FORMAT ( 1 H 1 ) 

16 FORMAT (54H THE GUESSED AMPLITUDE OF RESPONSE « , 

1 E14. 7 ,/ ) 

17 FORMAT ( 5(E14.7*6X) */) 

IS FORMAT ( 5X6HO ( L6 ) = * 1 7 X3HHU = * 1 2 X6HM ( SLUG ) = • 1 4X6HD ( 1 N ) = * 8X 1 2HW (RAD/ 
1 SEC)*) 

19 FORMAT (F10.b*7( E14. / *ix ) */ ) 

20 FORMAT (2X3HHA=*5X*6HFI(1> = *9X* 6HF I 13)=*9X3HFA=*12X*6HFI (7)=*9X*3H 
1DA=* 13X »4HDA0= ,11X*2HF= ) 

21 FORMAT (54H THE CHARACT E- 1 ST I C F REQUENCY 
1 E 14 . 4 * / ) 

22 FORMAT ( 54H THE NOND 1 MENS I ONAL I Z t U FREQUENCY 
1 E 1 4C4 * / D 

23 FO-MAT (+5,5X, ( 7F 10.5) ) 

24 FORMAT ( b4H iHt 2EKO— ORDtR ARiPl i T UOE OF RESPONSE * 


0 r IGIMa t „ 

***i«ss? 
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1 E 14 • 7 » / ) 

Zb FORMAT (54H THE NOND+MENSIONAl I ZED UPPER RESPONSE® » 

1 E 1 4 • 7 * / ) 

26 FORMAT I54H THE NGND I MENS ION AL I ZED LOWER RESPONSE* » 

1 E 1 4 • 7 * / ) 

BKS = C1> BCS*C2» AN1*N1=2«5 s AN2=N2*2*5, TOE IS THE TOLERANCE OF 
ERROR AALlOWED FOR THE SOLUTION, bMASS IS IN SLUG « AND DELTA IS 
STEP IN INCHES* 

AG ( I ) » AND AOG { I ) ARE I SERIES OF THE GUESSED VALUES OF AMPLITUDES# 
A AND AO RESPECTIVELY* ON ( 1 ) ThE GUESSED VALUES OF THE DYNAMIC 
LOAD IN LbF * WN(I) ARE SERIES OF THE NOND I MENS I ON AL I ZED FREQUENCY 
TO BE USED IN THE CALCULATION* 

-EAD 12 » BKS > BCS. ANl » AN2 
-EAD 12# HO 9 DM ASST DELTA! TOL 
READ 23. I G . < AG ( I ) * I - 1 * I G I 
-EAD 23. IG. I AUG ( I ) > 1 = 1 # +G ) 

READ 23. LQ. (QN( I) . 1 = 1 ,LQ> 

I A=0. 

DO 110 IQ=WLQ 

-EAD 2 3. LW*(V/N(I)tI=l*LWJ 

PRINT 15 

+AA=0C 

QS=QN( IQ ) 

AMASS® BMASS/ 12 <, 

AM DE L@ d e l T a?>$am A s s 

H0N1=H0**AN1 
&|ON2=HO*MAN2 
HOPNl*HONl*HO 
HOPN2~HON2*HO 
WS2*BKS*AN 1/AMDEL/HOPNl 
WS=SQRT 1 WS2 ) 

P-+NT 21. WS 
DO 110 +W31.LW 
+ A = + Q— 0 9 


WB=WN( I W) 

W r=WB*WS 
PRINT 10 

PRINT 17 * OS » HO s BMASS, DtuTA, W 

PRINT 22, W6 

Wb2=W6*WB 

WO-WS 2 * AMDEL 

Q=QS/WO 

HOMAX=HO*O.P9 

PI =3* 1416 

BK=BKS/WO 

BC=BCS*WB/WS/ AMASS 
PH01*PI*H0N1 
PH02-PI *HQN2 
3X*BK/PH01 
33=BC/PH02 

THE AMPLITUDE OF A for A0*0* 

IF (W.GT.u. ) GO TO 50 
I A A = 1 

AW0 = 0*P I *HON i / b K. 


ORIGINAL PAGE IS 
OF POOR QUALITY 


IS ESTIMATED BY LINEAR INTERPOSITION* 



IF ( AW0-F1A(50) ) 41, 42 * 43 

42 A~ 0 * 5#H0 
A0 = 0. 

GO TO 30 
41 11=1 

+ 2 = 30 
GO TO 44 

43 +1=30 
12=100 

44 13=1 11+12 )/2 
10 = 12 - 11 + 0*1 

+F (+D-1D 45* 45* 46 

46 IF ( AW0-F1A ( 1'2 ) ] 47* 49 * 48 

47 12=13 

GO TO 44 

48 11=13 

GO TO 44 

45 Al= + 1- 

A= ( A I + ( AWO-F1A (Il|)/<FlAU2)-FlA(+i)) >/100**HU 
A0=0. 

GO TO 50 
4 A = + 3/100* MHO 
A0 = 0 • 

50 CONTINUE 

GUESS THE VALUES OF A AND AO. 

IF ( IAA-1 ) 1*6,6 

1 CONTINUE 
+A=+ A+ 1 

IF ( IA.GT. IG ) GO TO 110 
A = AG( I A ) 

AO = AOG( I A) 

6 CONTINUE 
— +NT 16T A 
I T ER = 0 . 

103 CONTINUE 

HA05={ l.+AQ/HO>**2.5 
+A07= ( l.+AO/HO)**3*5*HO 
A+O3A/U+O+A0D 
HA=AH0*100. 

ha=ha+ 1 • 

I HA=HA 

+ A+ 1 

+ 1 O++AH+0Q A- + HA } * ( F 1 A l I HA 1 ) -F 1 A t I HA ) ) 

F + I2D3F2AU+HA ) + ( HA-+HA ) *UF2AU++A1 )-F2A( I HAD ) 

FI (3) =F3A( I HA I + < HA- 1 HA ) * ( F 3 A ( l HA 1 ) -F 3 A < I HA ) ) 

F 1 ( 4) =F4A { +HA > +U+A-1HA) * l F4A ( IHA1 >-F4Al I HA ) ) 
F+( 5) 3F 5A t + HA ) + { HA- + HA ) *UF5A( IHA1D-F5AI IHAU) 
FII6) =F6A( IHA) + ( HA-IHA) * ( F6A( I HAi )-r 6A( IHA ) > 

FI ( 7) =F7A ( I HA ) + ( HA- I HA ) * ( F7 A I I HA 1 ) -F7 A ( I HA ) ) 

FI (8)=F8A( IHA) + (HA-I HA ) * l F8 A l I HA 1 ) -F8 A ( 1 HA ) ) 
FI ( 9 ) =F 9 A ( I HA ) + ( HA- I HA ) * ( F9 A U H A 1 ) -F 9 A ( I HA ) ) 

F M 1 ) = FMD/HA05 
FI t 2 > = FI t 2) /HAG 7 
FI ( 3 ) = FI (3)/ I AO 5 
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F It 4 ) = FI ( 4 ) /HA07 
F I ( 5 ) =-F I ( 3 ) /HA07 
F I ( 6 ) =-F 1(6) /H AO 7 
FI ( 7 ) - FI ( 7 ) /H AO 5 
FI ( 8 ) = FI ( 8 > /HA07 
F I ( 9 ) *-FI ( 9 ) /HA07 

F = ( FI ( 1 ) * B 1-A* WB2 ) * *2 + ( F+ ( 3 ) *B3*A ) **2-G*Q 
F A 1 = ( F I ( 1 )*bl-A#WB2 ) * ( FI t 2 ) *bl-‘*lb2 ) *2* 

FA2= ( F I ( 3 J*B3*A) * ( F I ( 3) *B3 + F I ( A )*B3*A) *2« 

F A = FA 1 +FA2 
G = F I ( 7)-2.»PI 

FA0 = 2.*(B1*FI < l)-A*WB2)#ai*FI ( 5 ) + ( B3*A ) **2*F I < 3 ) *F II 6 ) * 2 » 

DEL-FA-m-F I ( 9 ) -F AO*F 1(8) 

0A= ( -F*F+ ( ) +G*FAO ) /DEl 

OAO= ( F*F I ( 8 ) -G*FA ) /DEL 

IF ( ABS(DA) .GT.CU5) GO TO 1 

AO” AO+DAO 

A= A+DA 

P-+NT 20 

PRINT 19 9 HA » FI(l)* FI ( 3 ) * FA* FI (7)* DA* DAO * F 

PRINT 10* FAO » DEL* G 

I TER= I T ER+ 1 

IF ( A.Lt.U. ) GO TO 1 

A0B-AB5 ( AO ) 

AA0=AB5( A-AOB) 

IF ( A AO * GT • HOM AX ) GO TO-1 

IF ( ITER.GT.l'; ) GO TO 1 

IF (ABS(DA) .G7.T0L) GO TO 103 

IF < ABSIDAO) •C.T.TOL) GO TO 103 

R=FI ( 3)*BC*A/PH02/(F1 ( 1 )*BK/PH01~A*WB2) 

ALPHA=ATAN ( R ) 

ALPHA=AuPHA*160./3. 1416 

PRINT 24, AO 

PRINT 13, A 

AU=A+AO 

ADOWN=A~AO 

PRINT 25, AU 

PRINT 26, A DOWN 

PRINT 14, ALPHA 

+aa=+aa+i 

110 CONTINUE 
END 


saasss? 
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Card 1 Format (8F10.5) 

BKS - Value of c., for the stiffness of the gas film force 

In lb f . 1 

BCS - Value of c« for the damping of the gas film force in 

lb f /ips. 

AN1 - Power n 1 for the stiffness force in terms of gas film 

thickness. 

AN2 - Power n ? for the damping force in terms of gas film 

thickness. 

Card 2 Format (8F10.5) 

HO - Normalized gas film thickness at equilibrium. 

BMASS Mass of the step ring in response In slug. 

DELTA Step depth of the pad in inches. 

TOL Convergence factor for the amplitude iteration. 

Card 3 format (15, 5X, (7F10.5)) 

1G - Total number of the guessed amplitude A of the 

response. 

AG(X) - Array of the guessed amplitude A of the response. 

Card 4 Format (15, 5X, (7F10.5)) 

1G - Total nuidber of the guessed amplitude Ao of the response. 

ADG(I) - Array of the guessed amplitudes Ao of the response. 

Card 5 Format (15, 5X, (7F10.5) 

LQ - Total number of the force excitation to be investigated 

in the production run. 

QN(I) - Array of amplitude of force excitation in lb f . 

Card 6 Format: (15, 5X, (7F10.5)) 

LW - Total number of the normalized forcing frequencies r.o 

be investigated. 


WN(I) 


Array of normalized forcing frequencies 
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If the sixth input statement Is located within the do loop of 
IQ «r l f L.Q, LQ sets of card 6 are required. 



I 
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PROGRAM RbRK.IT ( INPUT »OUTPUl , PUNCH * T AP E9 9 ) 

COMMON Y * OY » AT AbL * RTABL* 1FVD* X* DX*W* bX* bC* Q* AMASS* DELTA 
COMMON HO* TlAST* AY* FAC* Wb2 * HON1* HON2 * HQPNl* HOPN2 * AN 1 * AN2 
COMMON XA, TA* XT A * IA, IA2* OX l * WLb* Wd* oLOPE* TLW * wbO* XO 
DIMENSION Y ( 2 > * DY(2)» A TAbL ( 2 ) * -TAoLC2)* WO-JU18) 

DIMENSION X A ( 1 000 ) » TA(lCOQ)* XTA(IOOO) 

EXTERNAL, deriv, CNTRL* graph 

11 FORMAT ( 8F 1 0 • & ) 

12 FORMAT ( 5 AH MASS OF THE StAL (SLUG) * 

1 El4*3,/> 

13 FORMAT ( 5 AH DEEPTH OF STEP 0 1 SCUNT I NO I T Y (IN) * 

1 E14.3,/} 

14 FORMAT ( 54H NOND I MENS I ON Al I ZED EUUILlbKIUM POSITION* HO * 

1 E 14* 3 » / ) 

15 FORMAT { 54H FREQUENCY OF PlRIOuIl FORCE* W (RAD/SEC) » 

1 E14.3*/) 

16 FORMAT (54H AMPL. I TUDE OF THE PERIODIC FOKCc APPLIED ILLS) * 

1 E14. 3,/ ) 

17 FORMAT { S4H NGN 0 1 MENS lONALlZcD DAMPING CucFF IC itNT OR gAS F I LM 0 * 
1 E14.3,/) 

18 FORMAT (6uH NOND 1 MENS 1 ON Al I ZED MULTIPLE OF STIFFNESS * UlSPLACEME 
1NT * K*E14#3»/> 

19 FORMAT ( 54H NOND I MENS 1 ON Al 1 Z t D D i SPLACEMcNT » 

1 E 14* 3 * / } 

20 FORMAT (54H THE MULTIPLE Of TIME ANU FREOUtNCY (RAD) * 

1 E 14* 3 » / ) 

21 FORMAT ( 3X7HDfc.GREE**4X, 1 lHPHASt ( RAD ) = , 7X 1 3HD I SPL ACEMENT = » 1 1 X9HVEL0 
1C I TY= * / ) 

31 FORMAT ( 7X » 2 3HDATA ( XA ( I A J * I A* 1 * 70 ) / * 2 I L l A, 6 • 1H * ) ) 

32 FORMAT ( 5X * 1H2 , IX *4 ( E14, 6 * 1H » ) ) 

33 FORMAT (7X*22HDATA ( T A ( I A ) » I A= 1 , 9L ) / * 4 ( F 9 • a * 1H • ) ) 

34 FORMAT l 5X»1H3 *1X*6 l F9. 5* 1H* )) 

36 FORMAT ( 7X26HUA T A ( XT A ( I AP ) * I AP = 1 * 360 )/*2(El4*6*lH*n 
36 FORMAT ( 5X I H2 * 1 X * 4 { E 1 4. 6 * 1H * ) ) 

40 FORMAT (S2H ERROR RETURN* DX = Q */) 

41 FORMAT (52H NORMAL RETURN */) 

42 FORMAT (52H ERROR RETURN* VAKlAbLt INTERVAL MODE ONLY , / • 

TLAST IS THE FINAL TIME OF THE INTERVAL IN INTEGRATION* 

T LW IS The final time that The gMOOTH variaiion OF FREQUENCY END g • 
WLB IS The FINAL NONDIMENSIONAl IZED FREQUENCY TO BE JSED IN THE 
TECHNIC OF SMOOTH VARYING FREQUENCY. 

Y t 1 ) * Y 1 2 > » X ARE INITIAL CONDITIONS OF X, DX/DT, T RESPECTIVELY* 
BKS = C1. BCS=C2. DELTA + S STEP IN INCHES. OS IS DYNAMIC LOAD IN 
Lb* AM IS MASS Or THE PAD IN SLUG* Wb lb THE NOND IMtNS I UNAL I Z tD 
FREQUENCY* AN 1 = N 1 = 2*6* AN 2 = N2-Z* 5 0X1 Ig Tnc TIME INCRtMtN 7 U 
READ 11* T L A S T * T LW * WL3 
READ 11* Y ( 1 ) » r ( 2 ) » X 

READ 11 , BKS * bCS * QS, AM. DELTA* HO, WB 

0+ 

H0N1=H0**AN1 
HON2=HO**AN2 
HOPN1=HON1*HO 
HOPN2=HON2*HU ■ 


OBIGINAL PAGE IS 
OF POOH QUALITY 
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AMAS$*AM/12. 

AMDE AMAS S*UEl T A 
WS2-BKS *AN1/ AMDEL /H0PN1 
W3=SUR T ( WS2 ) 

W =£ WB*W5 
Wb2=WB*WB 
WL=WLB*WS 
xo= X 

SLOPE* (WLB-WB ) /WB/( TLW-XO) 

WB03WB 

WO=WS2*AMDEL 
Q*QS/WO 
8K=BKS/WO 
BC=»BCS/WS/AMASS 
PRINT 12, AM 
PRINT 13, DELTA 
PRINT 14, HO 
PRINT 13 , W 
PRINT 16, 06 
PRINT 17, BC 
PRINT 16, BK 

DEFINE THE INITIAL VALUE OF Y(I) 

XA<X>— Y(l) 

TA( 1 )=X 
X T A ( 1 ) = Y< 2 ) 

DEFINE X AS Y(l), AND DX/DT AS Y(2) 

the following is prepared for runoe-rutta numerical integration, 

P I *3 • 141692633b 
DX*DX I *P I / 13 0 • 

FAC=10.**U.2 

ntry=i 

N = 2 

I F VD~ 1 
I BKP* 1 

ATABL ( 1 )~0.00l 
AT ABL ( 2 ) = 0 • 0 0 1 
R TABL ( 1 )=0.001 
R T ABL ( 2 ) = 0 • 0 0 J. 

1 A2 = 1 
I A= 1 

PRINT 21 

CALL RKS3 (DERI V,CNTKL» Y,DYt AT AbL,RlAbL, WORK, X,DX,N, I FVO» I BRP »N i ,v 
1 I ERR ) 

IF (IERR-u) 4 r b» 6 
4 CONTINUE 
PRINT 40 

GO TO 3 
b CONTINUE 

PRINT 41 * 

GO TO 3 
6 CONTINUE 
PRINT 42 
3 CONTINUE 

THE Ful*- 0W l ,\vi ! .1 FOR The Pi_OfS OF X AND DX/DT V.S# T* 
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CALI NAMPlT 

call 3 C A L £ (XA*5.U*1A*1) 
call sc all ( ta.ao.o* i a, i ) 

CALL STALL <XTA.3.0*IA*1) 

CALL AXlb(O.O»O.0WHT V A LUt * 7 * 40 . 0 * 0 . 0 * T A <1 A+ 1 ) * T A (I A + 2 ) ) 

CALL AXIS < 0 . 0 , 0 . u*7HX VALUE *-7 *8. *90.0 *XA ( I A4 1 } *XA( 1 A + 2 ) > 

CALL AXIS ( 1 • C * 0 • 0 * 8HXT VAlUL * -8 * 8 • 0 * 90 • 0 » XT A II A+ 1 ) * XT A < I A+2 ) ) 
CALL LINL (TA,XA*1A, 1*1*0) 

CALL LINL I TA,X I A * I At 1 , 1 ,3 > 

CALL SYMBOL ( 3 .0 t9.0 tU.20* 16HPL0T OF X V.S. T,Q.Q,i6) 

CALL SYMBOLIC. 0* 1 0 . 0 *0 . 2 0 * 20HPLOT OF DX/DT V.s. T*0.Q*15) 

CALL ENDPlT 

THE FOLLOWING IS FOR The PHASE PLOT OF X V.S. DX/DT. 

R I A- I A 
WIDTH^IU. 

WSPACE= 10. 

CALL E I UPC (0 * >0* , 0 . ♦ 1 . *—20. * 2 . *WIDTH*WSPACE) 

CALL P AK AM ( GRA PH * 1 . »1.,RIA) 

CALL enusurf 
END 


UB-OUTINE PER IV 

DIMENSION Y<2)* D Y(2>* ATABL ( 2 ) , -TABL(2)t W0RM18) 

DIMENSION X A ( I 0 0 U ) * TA(IOOO)* XTA(IOOO) 

COMMON Y, DY * A1ABL* RTABL* 1PV0* X* DX.W* BK* BC* Q, AMASS* DELTA 
COMMON HO* T LAST , AY* FAC* WB2* HON I * H0N2 , H0PN1* HORN 2 » AN 1 * AN2 
COMMON XA> TA» XTA* IA* IA2* DX I * WLB* to/B* SLOPE* ILw, WBO* XO 
DY { ]. ) = D Y ( 1 ) / Dl * DY ( 2 ) =DY ( 2 ) /DT. 

DY ( 1 ) =Y < 2 ) 

DY < 2 ) = ( BK/HON 1 +(J*CuS ( X ) - BR/ ( HO-Y l 1) ) **AN 1-BC*Y ( 2 J *WB/ ( HO- Y ( 1 ) ) **AN 
12J/WB2 
RETURN 
END 


SUBROUTINE CNIRL(NTRY) 

DIMENSION Y(2:, DY ( 2 ) * ATABL12)* RTABL(2) 

DIMENSION XA(TOOO)* 1 A ( 1000) * XTA(IOOO) 

COMMON Y* DY * ATABL* RTABL* IFVD* X* DX*w* OK.* BC» U* AMASS* utL T A 
COMMON HO* T LAS T * AY* FAC* WB2* HONi* huN2 , H0PN1* H0PN2 » AN 1 * AN <• 

C 0 MMON XA* TA » XTA* i A * I A Z * DXi* WLG* rtt3* _ L’uP E * T L W * WBO » 

31 FORMAT ( SX * I 3*4 ( E 14.6 *6X > > 

34 FORMAT (34H THE D 1 SPL AC EMt NT IS OUT OF RANGc ) 

if y ( i > is less than -4.* or vm is greater than ho* terminaie 

THE NUMERICAL INTEGRATION. 

AY=ABS( Y( 1 ) ) 

IF(AY.GT.4.U ( G0T02 
. If { YU ) --H0 ) 1 , 2 , 2 

. CONTINUE 
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PRINT 5 A 
NTRY=2 
1 CONTINUE 

IF ( X-TlAST ) 3 , 3 » 4 

A NTRY- 2 
3 CONTINUE 
I A 2 = IA2 + 1 
I ADX= I A2*DX I 
I A= I A2/A+0.00C1 
X T A < I A ) =Y ( 2 ) 

X A ( I A ) ~ Y t 1 J 
TAUA1A 

PRINT 51 * IADX* X* Y ( 1 j » Y(2), WB 
II (Wb-WLB) 11 , 12* 12 

11 WB=WBO* ( 1 .+SLOPE* (X-XO J ) 

GO TO 13 

12 WB-WLB 

13 CONTINUE 
WB2=WB*WB 
RETURN 
END 


SUBROUTINE GRAPH l T * XP * YP * ZP ) 

COMMON Y, DY, ATABL, RTABL* I FVD* X* DX*W* BK. * 8C* Q* AMASS* UElTa 
COMMON HO, T LAST » AY* FAC, Wb2 , HON1 * HON2* HORN 1 * HOPN2 * AN 1 * A NZ 
COMMON XA* TA, XTA, I A, I A2 * DX I , WLB, Wb* SLOPE, TLW * WBO, Xu 
DIMENSION Y ( 2 ) , DY { 2 ) * ATA6L12)* R T A b L ( 2 ) * WORKI1S] 

DIMENSION X A { 1 00 0 ) , TA(IOOO), XTA(IOOO) 

NP-T+ • 0 0 0 1 
XP=X A(NP)*5. 

YP -U • 0 

ZP=XTA ( NP ) *5, 

RETURN 

END 
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Card 1 Format (8F1Q.5) 



TLAST 

- 

Final normalized time in radian for the Integration. 


TLW 

- 

Final normalized time In radian for the changing of 
V during the Integration. 


WLB 

- 

Final normalized forcing frequency. 

Card 

_2 

Format (8F10.5) 


Y(l) 

- 

The Initial value of X. 


T<2) 

- 

The initial value of X, i #e . dX/dT. 


X 

- 

The initial time, for the integration. 

Card 3 

Format (8F10.5) 


BKS 

- 

Value of c. In lb_. 

x r 


BCS 

-■ 

Value of C 2 In lb^/ips. 


QS 

-• 

Value of q in lb.* 
t 


AM 


Mass of the ring in response in slug. 


DELTA 


Step depth, 6 in inches. 


HO 

- 

Normalized equilibrium gas film thickness. 


WB 

- 

The initial normalized excltational frequency. 

Card 

£ 

Format (8F10.f) 


AMI 

- 

Value of n, which Is 2.5 in the case being investigated 
here. 


AM2 

- 

Value of n« which is 2.5 in the case being investigated 
here. 


DXI 

- 

Time increments during the fixed Interval of the numerical 


Integration. 
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PROGRAM RSTABl INPUT * OUTPUT > 

0 I MEN S l ON SA(30> * SB(30), SC(3U), SD(50)» 5E13C)* SF(3Q) 

DIMENSION PO ( 30 ) * PCR130}* PC I I 30 ) * DX<30), C(IOO), AI70.70) 
DIMENSION C II 2 * 3 * 30 ) » SUBV l 30 ) , V I < 20 ) 

DIMENSION XXI 100 ) , B ( SO , 3 ) * BH( 30 ) * AH (30)* CH( 50 1 

COMMON Bl* ALAM* HI* TOL* NL * N2* B2 * NLl* N21* NP * NP1* IML2* Nl3 
COMMON H2* H 1 3 * H23, HSUM . 

COMMON SA> SB* SC* 5D* SE* SE 

1 1 FORMAT (4F10. 4*3151 

12 FORMAT (8F10.3) 

13 FORMAT (54H THE THRESHOLD FREQUENCY OF THE STEP SEAL 
1 E 14 • 4 * / ) 

14 FORMAT (54H THE SQUEEZE F I lM PARAMETER 
1 E 1 4 • 4 * / ) 

13 FORMAT (34H THE BEARING NUMBER 
1 E 14 • 4 * / ) 

16 FORMAT (34H MASS 
1 E 14. 4 */// ) 

17 FORMAT ( 54H RATIO OF MIN CLEARANCE HEIGHT TO D IFF. IN CLEARANCES 


1 E 1 4 • 4 , / ) 

18 FORMAT ( 5 4H THE 1ST ORDER PRESSURISED LOAD * 

1 E14.4,/) 

19 FORMAT 1 l 8E14. 4) , / ) 

20 FORMAT ( 34H THE POINTS ARE TOO FEW ) 

21 FORMAT ( 54H A AND/OR B IS OUT OF RANGE OF TABLE ) 

22 FORMAT (34H THE GUESSED THRESHOLD FREQUENCY OF THE STEP SEAL * 

1 E 1 4 . 4 * / ) 

23 FORMAT (52H THE GRID DIFFERENCE */) 

24 FORMAT (52H THE ZERO ORDER PRESSURE DISTRIBUTION */) 

23 FORMAT (34H THE VALUE OF SMALLNESS * 

1 E 1 4 • 4 * / ) 


26 FORMAT ( 54H LOAD DUE TO THE IMAGINARY PART UF COMPLEX PRESSURE . 

1 E 1 4 * 4 » / ) 

27 FORMAT ( 3X * I 3 * ( 7F 10 . 4) ) 

29 FORMAT ( / , I 8 El 4. 4 , / ) ) 

30 FORMAT (32H THE IMAGINARY PART OF COMPLEX PRESSURE DISTRIBUTION,/} 

31 FORMAT ( 5 2 H THE REAL PART OF COMPLEX PRESSURE DISTRIBUTION */) 

32 FORMAT (34H THE RATIO OF WIDTH OF STEP WITH HEIGHT Hi , 

1 E 1 4 * 4 * / ) 

READ 11* Bl» AlAM* HI* TOL, NL, N2* NV 
READ 12 * ( V I ( 1 ) » I = 1 * NV ) 

B2=1.-B1 
NL 1 - NL — 1 
NL2=Nl~2 
NL3-NL-3 
ND^NLZ+NLZ 
N2 1 "N2- 1 
NP = N2+ 1 
NP1=NP+1 
H12=H1**2 
H13 = H1**3 ( 

H2 =H 1+ 1 * 

H2 2'=H2 * *2 
H23~H2**3 
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READ 12. <DXU I .1=1, NL1) 

READ 12 » (POC I ) , 1=1, NL) 

HSUM=H13+3.*H1*H2*(H1+H2 )+H23 
CALL NEWR(PO.DX) 

CALL INTEGt 0..1- tDX.PO , Nl » WLQAD , C I , I ERR ) 

PRINT 32, 61 

PRINT 15, ALAM 

PRINT 17, HI 
PRINT 14, SEG 
PRINT 23 

PRINT 19, ( DX ( I ) . I * 1 .ML! ) 

PRINT 24 

P R I N T 19 , (POM), 1 = 1 . NL ) " “ 

PRINT 25, TOL 
Al=-ALAM#Hl/2. 

A2=~ALAM* (H1+H2 ) /4. 

A3=-ALAM*U2/2. 

THE COEFF I C i Lf'lT b OF SMALL A, d, C, 0, E, F AT EACH GRID POINT akl 
CALCULATED IN THE FOLLOWING AS SA , ... ETC. 

DO 2 1=2, N21 
SA( I ) =H13/DX( 1 ) 

SC ( I ) =H 13/DX{ J-l ) 

SB ( I ) =-SA ( I ) -SC ( I ) 

SD1= ( PO ( I + 1 ) **2/DX t I ) -P0( I l./DX ( I J+l./DXt 1-1) )+P0( I-l)**2 

1 / DX ( 1-1) ) #3 • +■ HI 2/4« 

SD2=ALAM* (PO ( i +1 )-P0 H-ill/2. 

SDt I) =SD1-$D2 

SE( I ) =H 1* ( DX ( 1 )+DX( I — X ) ) /2* 

SF(I ) =( DX( I ) +DX ( 1-1) ) /2m 
AH(I)=A1 

chi d=-ai 

BH( I ) =0. 

2 CONTINUE 

SA(N2 ) =HSUM/8.,/DX (N2 ) 

$C(N2)=H13/DX i! N 2 1 ) 

SB(N2)=-SA(N2)-SC(N2) 

SD1 = ( ( PO ( NP ) **2-PO(N2 )**2 ) * ( H12+H 1*H2+H22 ) /DX(N2)+(PQ<N21)**2- 
1 PQ(N2)**2)*3.*H12/DX(N2l)>/4. 

$D2=ALAM*(PO(NP)-PO(N2l) ) / 2 • 

SD( N2 ) =SD1-SD2 

SE ( N2 ) = ( (H1+H2)*DX(N2)+2.*H1*DX(N21> )/4. 

SF(N2) =(DX(N2)+DX(N21> ) /2. 

AH(N2)=A2 
CH(N2 ) =-Al 

BH(N2)=AH(N2)+CH(N2) 

DO 3 I =NP 1 ,NL.L 
SA( I )=H23/DX( [ ) 

SC ( I ) =H23/DX { t - 1 ) 

SB ( I ) =-SA ( I ) -SC ( I ) 

SDl-( P0( 1 + 1 ) * *2/UX< I }-P0< I >**2*< l./DX ( I ) + l./l>X( 1-1 ) )+P0( 1-1 )**2 
1 / DX ( I - 1 ) ) *3 • *H 2 2 / 4 * 

SD2=ALAM*(P0< 1 + I )-P0( I-i ) ) /2. 

SD< I ) =SD 1-SD2 

SE (I ) =H2*(DX i i ) +DX U -1 ) ) / 2« 



- 120 - 


SF( I) = (DX(n+DX< I~l> ) /2m 
AH ( I ) = A 3 
CH(I)=-A3 
8H( I ) = 0 • 

3 CONTINUE 

SA ( NP ) -H23/OX ( NP ) 

SC ( NP ) =H$UM/8*/DXt N2 ) 

SB( NP ) --SA t NP ) -SC t NP ) 

SD1=( (P0(NP1 )**2~PO(NP)**2)*3.*H22/UX(NP)+<PO{N2)**2-PO<NP>**2)* 
1 (H12+Hl*H2 +H 2 2 J /OX ( N2 ) ) /A# 

SD2*ALAM*(P0(NP1 )-P0tN2> >/2. 

SDt NP J =SD1-5D2 

SE(NP>*( (H1 + H2)*DX(N2H-'2.#H2*DX<NP) ) /4. 

SF ( NP ) = ( DX ( NP J f DX ( N2 ) )/2. 

AH ( NP ) - A3 
CH ( NP ) --A2 

BH(NP)=AH(NP)+CH(NP) 

CALCULATE ELEMENTS OF (A) AND (C) IN THE EQUATION AX*C* 

B MATRIX IS USED TO SAVE A(I*J) WHICH ARE INDEPENDENT UN THF 
DO A J = 1 , ND 
DO A 1=1, ND 
A ( I * J ) = D. . 

4 CONTINUE 
V=VI ( 1 ) 

PRINT 22* V 
AMASS = Wl.OAD/ V#*2 
PRINT 16* AMASS 
DO 5 J = 2, NL3 

I - J + l 
JN=NL2+J 
t N = JN+ 1 

A( J* I ) =SA( I )*P0( 1 + 1 ) +AH ( I ) 

A ( JN * I N ) = A ( J * I J 
A ( J * J ) =SB ( I ) *P0{ I ) + BH ( I ) 

A( JN*JN)=A( J.Jj 
A(J*J-lJ=SC(n *P0 t 1-1 >+CH ( I ) 

A ( JN * JN- 1 } =A ( J , J-l ) 

A( J* JN } =SE ( I ) *V 
A ( JN* J ) =- A t J » JN > 

C( J)=-SDt I ) 

C( JN ) =SF ( I ) #P0( l ) *v 
B(J*3)=A(J,J+1) 

B< J,2) =A{ J, J) 

B I J * 1 I -At J * J-l j 

5 CONTINUE 

A ( 1 * 1 ) =SB ( 2 ) *P0( 2 )+BH( 2 ) 

A ( NL 1 , Nul )=A( 1,1) 

A ( 1 * 2 ) = SA ( 2 ) *P0 < 3 ) +- AH ( 2 ) 

A ( NL 1 »Nl. ) = A ( 1 , 2 ) 

At 1 *NL1 ) -SE t 2 > *V 
A(NL1*1 ) = -At 1 ,NL1 ) 

C ( 1 ) -“SD ( 2 ) 

C t NL1 ) =SF t 2 ) *P0( 2 ) *v 
A(NL2,NL_3)=SCtNLl)*PO(NL2) +CH ( NL 1 ) 



A ( ND » ND- 1 )-A(Nl2,NL3 ) 

A<NL2,Nl. 2)=SB(NL1 )*P0( NLl ) +BH t NLl ) 

A ( ND .ND ) = A ( NL2 , NL2 ) 

A(NL2,ND) =SE(Nul )*V 
A ( ND*NL2 ) =-A ( Nu2 *ND } 

C (NL2 ) =-SD { NLl ) 

C<ND)=SF(NL1 )*PO(NLl )*V 
B ( 1 * 1 ) =A( 1 , 1 ) 

B( 1 *2 ) *A( 1 t2 i 
B ( NL2 , 1 ) =A( NL2 ,NL3 ) 

B(NL2,2 )*A(NL2 ,NL2 ) 

DO 61 J s 1 »ND 
XX(J)=C(J) 

61 CONT I NUt. 

CALL DETECH A.XX.ND.DET ) 

IF (DET-U.) 1 u 8 i 109. 108 
109 5T0P 
108 CONTINUE 
PCI(1)=0. 

PCI (NL) =0. 

PCR ( 1 ) - 0 # 

PCR( NL ) *0, 

DO 7 J = NL 1 . ND 
PCIU-Nl3)=XX(J) 

7 CONTINUE 
PRINT 30 

PRINT 29. ( PC 1 ( I ) * 1 * 1 «NL ) 

CALL I NT EG ( 0. . 1. .DX.PCI .Nl.PGMI »CI . I ERR ) 
IF C I ERR — 1 1 113* 1 1 A » 116 
1 1 A PRINT 20 
STOP 

115 PRINT 21 
STOP 

113 CONTINUE 

PRINT 26, PSM1 
DO 69 J=l» NL2 
69 PCR ( J+ 1 } = XX ( J ) 

PRINT 29. ( PCR ( I ) . I = 1 * NL ) 

PRINT 31 

CALL INTEGtO. *1. # DX.PCR.NL tPRbM.CI.IERR) 

PRINT 18. PRSM 

AMASS*PRSM/V**2 

PRINT 16. AMASS 

DO 121 I V = 2 » NV 

V^VI (IV) 

PRINT 22, V 
DO 51 J = 1 »ND 
DO 51 1*1. ND 
AfUJhO, 

51 CONTINUE 

DO 52 J = 2,IML3 
I = J + 1 
JN-NL2 + J 
IN=JN+1 
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A( J*l) 

A(J*J) 3= y(J»2j 
At J*J+1 ) S B( J*3 ) 

At JN*IN)=At J*I ) 

A ( JN * JN ) =A ( J * J ) 

At JN* JN-l)=AtJ*J-ll 
52 CONTINUE 

At 1* 1 ) *BC 1 . 1 ) 

At 1 > 2 ) t 1*2) 

A ( NL 2 * Ni_3 ) =B < Nl2 * 1 ) 

A ( NL2 * Nu2 )=B(Nl2*2) 

A t NL 1 * Nu 1 ) = A ( 1 * 1 ) 

A ( NL1 *Nu ) S A ( 1*2) 

A(ND*ND~1 > =A t |Vl2*NL3 ) 

AIND*ND)-A(NL2*NL2) 

CALCULATE THE elements OF A AND C WHICH VARY WITH RESPECT TO V. 

DO 6 J=1*NL2 

I = J+l 

JN=J+NL2 

At J.JN)=SE( I )*V 

A ( JN * J ) = -A ( J * JN ) 

C( JN) =SF( I ) *PO( I )*V 
6 CONTINUE 

DO 62 J = 1 *ND 
XX(J)=CU) 

62 CONTINUE 

CALL DETEQt A*XX*ND*DET) 

IF IDET-O.) Ill* 112* 111 
112 STOP 
111 CONTINUE 

DO 8 J=NL1.ND 
PC I t J-Nu3 ) S XX i J ) 

8 CONTINUE 
PRINT 30 

PRINT 29* ( PC I { I ) » I * 1 *NL ) 

CALL 1NTEG( 0. . 1. »DX.PCI , NL * PSM2 , C I * I ERR ) 

IF ( IERR-1) 116*117. 118 

117 PRINT 20 
STOP 

118 PRINT 21 
STOP 

116 CONTINUE 

PRINT 26* PSM2 

DO 9 J = l, NL2 

9 PCR ( J+ 1 ) - XX t J ) 

PRINT 31 

PRINT 29, l PCR ( I ) * I = 1 *NL ) 

CALL INTEG( 0. ,1. *DX,PCR,NL,PRSM,CI ,IERR) 

PRINT 18, PRSM 
AMASS-PRSM/V**2 
PRINT 16* AMASS 
121 CONTINUE 
END 
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* * 


# 

* 


* 




If 

tf 


•ft 


* 

* 


* 


* 


SUBROUTINE I NTE6 { A »6 »H , F *NP ♦ VALUE *C * I ERR ) 

******#**#* # *** # #* <t## , <######## 

SUBROUTINE INTEG 


INTEGRATES THE NON EQUIDISTANTlY tabulated FUNCTION FiXUV) 
BETWEEN THE LIMITS A AND 6. 

A modified method of overlapping parabolas is employed. 

A SECOND ENTRY POINT i lNTEG2' 1^ PROVIDED FOR MORE THAN ONE 
INTEGRATION ON THE SAME DIVISIONS OF X, THIS SAVES THE TIME 
OF CALCULATING THE WEIGHTING FUNCTIONS, 

ARGUMENTS - 

A LOWER LIMIT OF INTEGRATION. 

B UPPER LIMIT OF INTEGRATION, 

X ARRAY OF ARGUMENT VALUES, MUST BE MONOTON I CALL Y 

INCREASING AND MUST bt DIMtNSIONED NP, 

F ARRAY OF FUNCTION VALUES, MUST BE DIMENSIONED NP. 

N p NUMBER OF POINTS. NP MUST bE GRtATtR THAN 3. 

VALUE RESULTANT VALUE OF THE INTEGRATION, 

C WEIGHTING FUNCTION PASSED TO THt MAIN PROGRAM 

FOR STORAGE. 

I ERR RESULTANT ERROR PARAMtTtR, 


* * * INTEG002 
♦INTEG003 

* I NTEG0Q4 
♦INTEG0Q5 
* INT EG006 
*INTEG007 
♦INTEG0Q8 
* INTEGQ09 
♦INTcGOlO 

* I NT EGO 1 1 
* INTEG012 
* INTEG0I3 
*INTEGQ14 
* INT EG01 5 

* I NT EGO 1 6 
* INT EGO 1 7 
* INTEG01S 
* INTEGQ19 
* INT EG020 

* I N TEG02 1 

* 1NTc6022 
*1NTEGQ23 

* I N T EGO 2 A 
♦INTEG025 
* INT EGO 26' 


* 


■* 


# 


■ * 

■* 


it 

* * * 


* 

■y 

# 

t 

y 


REQUIRED SUBPROGRAMS - NONE 


COMMON STORAGE - 

THE WEIGHTING FUNCTION C IS STORcD IN THE MAIN PROGRAM AND 
REQUIRES THE FOLLOWING DIMENSION STaTEMlNT WHlRE d.GE.NP. 
DIMENSION C t 2 * 3 * D ) 


ERROR INDICATIONS - 

IERR = 0 INDICATES NO ERROR. 

IERR = 1 INDICATES NP IS i_ESS ThAN 4. 

IERR = 2 INDICATES THE LIMITS OF INTEGRATION ARE OUT OF 
THE RANGE OF THE TAsLE. 


EDWARD G. TRACHMAN 8 JULY 1970 M.E. DEPT. 492-3640 

**** > *****#****# ##4( . ###### . lt<###<## 
DIMENSION X ( 30 : . F(30)* C(2,3 *d0) 

DIMENSION H ( 30 } * SUBV ( 5 0 ) 


NP MUST BE GREATER THAN 3 
IF ( NP * u.E * 3 ) GO TO 96 


CALCULATION OF INTERVALS OF X 

NH=NP-1 
X( 1 )=0, 


* 1 NT E602 7 

* I NT EG028 

* INT EG029 
* INTEG030 
* 1NTE6G3 1 

* IN T 1 6032 
*IN TEG033 
* INT EG034 
* INTEG035 

* IN I EGO Jt) 
* INTEG03 7 

* IN IEG038 
*INTEGQ39 
*INT E6Q40 
♦INTEG041 

* I NT E6042 
* INTEG043 

IN I EG045 
. i N i uvjOho 
INI toQ 4 *? 
I N f c. GD48 
I N I LG049 
iNTcGOSO 
I N T cGQbl 
INTEG032 
INT EGG33 



124 


DO 10 1*1 » NK 
10 X ( I + 1 ) =X ( I ) +H { I ) 

DO 20 1 = 1* NH 

IF ( I • EQ* 1 ) GO TO 15 

DEFINE COEFF I C I ENTS OF FIRST PARABOLA 

C( 1*1*1 I) )**3/(6.*HU-l)*(H( 1-U+HlI) .11 

C( 1 *2* I)=H( I )*( 3.*H( I~1 )+H l 1 ) ) / <6.*H(I-1 )) 

CI1*3*I ) = H t I ) * l 3 .*H < I - 1 ) +2 *#H ( I ) )/(6**(H(I-l)+H< I) ) > 
15 CONTINUE 

IF (I.EQ.NH) GO TO 20 

DEFINE COEFFICIENTS OF SECOND PARABOLA 

C ( 2 * 1 * I ) = H ( I )* ( 2 •*n( I )+3 **rl< i + 1 ) ) / ( 6.*lrt ( I > +HI 1 + 1 ) ) ) 
CI2*2»I)=H( I ) * ( h ( I ) + 3 • *H ( 1 + 1 ) ) / I 6 • *H ( i + 1 H 
Ct 2*3*1 )*= — (Ht I ) ) ** 3/ 1 6 • * H { 1 + 1 ) * ( H ( I ) +H { 1 + 1 ) H 
20 CONTINUE 

•ENTRY INTEG2 


INITIALIZE SUMMATION VARlAbLt 

value=o.o 

IF (B-A) 40*92,30 

B IS GREATER THAN A 

30 ALIM = A 
BLlM = B 
SIGN * 1*0 
GO TO 50 

A IS GREATER THAN d 

40 ALIM = B 
BLlM * A 
SIGN = -1.0 
50 NH=NP-1 

SETTING THE LOWER LIMIT OF INTEGRATION 

DO 63 1=1 »NH 
PARTA = 1.0 

IF (ALIM-X(I)) 61,69*63 
61 IF ( I • EQ* 1 ) GO TO 97 
ALIM = X( 1-1) 

PART A= ( X { I )-ALlM) / (X( 1)-X( I-l) > 

GO TO 69 
63 CONTINUE 
69 CONTINUE 


ORIGINAL PAGE IS 
OF POOR QUALITY 


lNTtG056 

INTEG057 

INTEGQ58 

INTEG059 

INTEG060 

INTEG061 

I NT EG062 

INTEG063 

INTEG064 

INTEG065 

I NTEG066 

INTEG067 

I NT C.G068 

I NT t G069 

INTEG070 

INTEG071 

INTEG072 

INTEG073 

1NTEG074 

1NTEG075 

INTEG076 

INTE6077 

INT EG078 

INTEG079 

INTE6080 

1NTEG081 

1NTEG082 

INT EGOS 3 

INTEG084 

INI LbQaS 

1NTEG0E6 

1NTEG087 

INTEG088 

INTEG089 

I NT t G090 

l NT cG09 1 

INTEG092 

INT LG093 

1 N V t G094 

1NTEG095 

INTEG096 

INTEG097 

INTEG098 

INTE6099 

1NTEG100 

TNTEG101 

iNT EG 1 02 

1MC.G1Q3 

I N T lG 1 04 

INT EG105 

INTEG106 

I N T r. G 1 0 7 
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* 

* 


* 

* 

* 


* 

# 

* 


* 

* 

* 


* 

* 

* 


SETTING THE UPPER LIMIT OF INTEGRATION 

00 73 1=1. NH 
PARTB =1.0 

IF (BLIM-Xl I ) ) 71.79*73 
71 IF ( I * EQ. 1 ) GO TO 97 

BLIM = X ( I ) 

PART8=(BLIM-Xf 1-1 M/(X( I )-X( M M 
GO TO 79 

73 IF ( I.EQ.NP) GO TO 97 

79 CONTINUE 

CALCULATION OF INTEGRAL OVER SU6 INTERVAL 

DO 80 1=1 » NH 
SUBV( I) s 0.0 

IF ( X ( I ) .EO. Al IM ) -3U8V1 I )=C (2.1*1 ) *F ( I )+C(2»2*l)*F(I+l)+C(2*3*I)* 
IF l 1+2) 

IF U’NH) 102* 103 . 103 
103 CONTINUE 

A0X = ABS < X < 1+1 i-BLIM) 

IF ( ADX.LT. l.E-5 , SUBV( I )=Ct 1*1*1 ) *F ( I - 1 ) + C { 1 * 2 * I) *F ( I ) +C ( 1 * 3 * I } * 
1 F ( 1 + 1 ) 

GO TO 101 
102 CONTINUE 

IF < X ( n.GT.ALlM.AND.XU+U.LT.OLlM) SUBV( i>= 0 . 5 *(C(l*ltl)*FtI-l) 
1+ ( C ( 1 * 2 * I ) +c ( 2 , 1 * I) ) *F ( n + ( C ( 1 * 3 . 1 ) +C ( 2 * 2 . I) ) *F i I + 1 ) +C ( 2 * 3 . I ) * 

2F( 1+2) ) 

101 CONTINUE 

IF (PARTA. NE.l. OR. PARTB. NE.l ) SUB V ( i )= PART A*P ART 6*SUbV ( I) 
CALCULATE THE FINAL VALUE OF T Ht INTEGRAL 

80 VALUE= VALUE+SUBV { I ) 

VALUERS IGN*VALUt 

SET ERROR PARAMETER FOR NORMAL RETURN 

92 I ERR = 0 
RETURN 

SEJ ERROR PARAMETER FOR TOO FEW POINTS 

96 I ERR = 1 
RETURN 



97 I ERR = 2 
RETURN 
END 


8; tfi-m T ER.-f.OR d&fctf -0+:' 1 AW* 0 


I NT EG 108 
I NT EG 109 
INTEG110 
INTEG111 
INTEG112 
INTEG113 
INTEGUA 
INTEG115 
INTEG11 6 
INTEG11? 
INTEG118 
INTEG119 
I NT EG 120 
1NTEG121 
INTEG122 
INT EG 123 
I NT EG 124 
INTEG123 


I NT EG 126 
INTEG127 


INTEG128 
INT EG129 
INT EG130 

INTEG13 1 
INTEG132 
INTEG133 
INTEG134 
INT £0135 
INTEG136 
i NT EG 1 3 7 
INTEG138 
INTEG139 
INTEG140 
INT EG 14 1 
INT EG 142 
INT EG143 
I NTEG144 

i n r eg 145 

l ;:T lG 146 
1 i. O i *+.7 

-iV ^ ! -i : "-i: 48 - 

1 'Fi rS (60 
‘ ' ■: -o i :> 1 
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SUBROUTINE DETEQ { A »B * N » DET ) 

D 1 MENS 1 ON A t 70 » 70 ) .B l 70 ) * 1PV0T < 70 ) 

DET-1.0 

DO 11 J-l.N 

II IPVOTJ J ) = 0 ♦ 0 
DO 121 I = 1 * N 
T = 0 • 0 

DO 51 J ~ 1 *N 

IF ( I P VOT t J ) — X > 21*51.21 
21 DO 50 IC S 1 »N 

IF ( IPVOT(K)-l) 31.50*50 
31 IF (ABStTJ-A3S(A( J.K) ) ) 41. 50 .50 
41 I ROW= J 
I COL~ 1C 
T = A ( J * K ) 

50 CONTINUE 

51 CONTINUE 

IF ( ABS ( T ) -1 • E-8 ) 131*131.55 

55 I PVOT ( IC0U=1 

IF { I ROW- I COL ) 61.61.61 
61 DET=-DET 

DO 71 L - 1 . N 
T = A ( IROW.L > 

A ( I ROW . t_ ) = A t ICOL.LJ 
71 AUCOL.O^T 
T -B ( IROW ) 

B( I ROW ) - B ( ICOi. ) 

B( ICOL ) =T 

81 TEMP=A( ICOL • ICOL ) 

DET=DET*TEMP 
A ( I COL . I COL ) = 1 • 

DO 91 L=1.N 

91 Al ICOL.u)=AI ICOL. L) /TEMP 
B l I COL ) S B l I COL ) / T EMP 
DO 121 uI«l.N 
IF t L 1 - I COL ) 101.121. 101 
101 T =A ( L 1 . I COL ) 

A ( L 1 » I COL ) «0 • 

DO 111 u=l.N 

III A(L1.L)=A<L1.L)-A(ICGL.LI*T 
B(L1)*B(L1)-B( I COL ) * T 

121 CONTINUE 
RETURN 
131 DET=0.0 
RETURN 
END 


SUBROUTINE NEWR(PO.DX) 

DIMENSION AA (50)* bbl50). CC150). DD(50>. EE15Q). FF t 50 ) * DX15GI 
DIMENSION A l 50 ) » Bt50). C<50). PO(50). DP150). F I t 50 > * FIDPI50.3) 
COMMON Bl. ALAM. Hi. TOL » NL » N2 » B2 * NL 1 ♦ N21. NP * NP 1 * NL2 . Nl3 



COMMON H2 * H13 * H23* HSUM 
COMMON AA, BB * CC* DO* EE* FF 

11 FORMAT (4F10.4*2I5) 

12 FORMAT ( 8F 1 0 • 5 ) 

13 FORMAT (6XE14*4,3{6X»El4.4) ,/ > 

1 A FORMAT { 17X3HFI=*3X17HDFI/OP ( J ) * J*l*2*3* / ) 

15 FORMAT (3X2HI=*2 2X3HDX* # 1 7X3HPO* * / ) 

16 FORMAT ( 5X15*2 (6XE14.4) ,/ ) 

17 FORMAT 154H PRESSURE 1$ SMALLER THAN Trlt AMBIENT PRESSURE 

18 FORMAT ( 5X15 ,3 (6XE14.4) ,/ ) 

19 FORMAT (1H1) 

DO 2 1*2* N21 

AA ( I ) *H23/2 • /DX ( I > 

CC( I) =H23/2 • /DX ( 1-1 ) 

BB (I ) ss- AA ( I >~CC< I ) 

FF ( I ) =Ai-AM*H2/2* 

DD( I) = ~FF i I ) 

EE t I ) =0 # 

2 CONTINUE 

AA ( N2 ) *HSUM/ DX ( N2 ) / 16 * 

CClN2)=H23/2./DXlN21 ) 

B6(N2)=-AA(N2)-CC(N2 i 
FF(N2)=FF(N21 ) 

D0IN2)*-ALAM*<H2+H1 >/4. 

EE ( N2 ) *DD ( N2 J + FF ( N'2 ) 

DO 3 I«NP1* NL1 
AA ( I ) *H 13/2 • /DX ( I j 
CC( I J*Hl3/2./OXf 1-1 ) 

BB ( I ) =-AA ( I >-CC< I) 

FF ( I )»AuAM*Hl/2. 

DO ( I ) *-FF< I ) 

EE( I ) *0. 

3 CONTINUE 

AA(NP)*H13/2*/DX(NPl 
CC(NP) =HSUM/DX ( N2 ) / 16 • 

Bb( NP) *-AA( NP >-CC( NP > 

DO(NP) «DD( NPU 

FF ( NP) =ALAM* (H1+H2 ) /4* 

EE(NP)=DD(NP)+FF(NP> 

104 CONTINUE 

DO 5 1=2* NL 1 
J = I-1 

A( J)=AA( I )*P0< I+lJ+DDTl 1 
B( J)*8B( I l*POl I ) +EE ( I ) 

CC J)=CC( I )*P0( I-1)+FF( I ) 

FKJ)*A(J)*P0(I + 1 >+B( J) *P0< I )+C ( J )*POl 1-1) 

FI ( J ) = -F I ( J) 

5 CONTINUE 

DO 6 I * 3 * NL2 
J=I-1 

F I DP ( J * 1 ) = CC( n»2**P0( 1-1 )+FF ( 1 ) 

F I DP ( J * 2 ) = 63 ( I )*2.*P0< I |+tU I ) 

F I DP ( J * 3 ) = AA ( I )*2**P0< 1+1 >+DU( 1 ) 

6 CONTINUE 
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hIDPU.2)=BB(2)#2**P0(2)+tL(2J 
p 1 UP ( X , j ) r A A ( 2 ) » i • * P 0 < 3 ) + D L) ( A ) 

F | DP ( NL2 > 1 ) = 0 0 l ML 1 ) *2 • l«L A > + h h t NL a ) 

F1DPINL2*2) - bo l Nl i ) *2 • *PO t NL i ) +ll l NL i ) 
CALL TULtO ( h I DP * F 1 » DP . Nl A ) 

DU 7 I =2. NLl 

PO ( I ) =pu ( I )'+DP ( 1 - 1 ) 

Ih (PO( 1 ) - 0. ) lol * 10 1. 7 

1 J 1 PRINT 17 

DU 23 I P= 1 * Nl.l 

PK 1 NT 18 * IPjuXI IF) »PU( IP) » up < 1 P -1 ) 

23 CONTI NUL 
STOP 

7 CONTI N'JL 
R M A X = 1 Ou/2 . 

DO 8 K = 1 » N E 2 
R = DP( K ) 



1 E i R-RMAX ) 8 % 

8 . I 'J 2 


102 

R M A X = R 



8 

CUNT I NUL 
I h { RMAX-IOL ) 

10 3* 1 U 3 

, 104 

103 

CUNT I NUt 
PRINT lb 
DU V I P = 1 *NL 1 
PR I NT 16 * IP, 

U X ( IP) . 

PU( IF) 


9 CUNT I NUL 
PRINT IV 
RETURN 
END 


SUHROUT1NL lULtu lA.O.X.N) 

S'JBRUU T i Nl oJlVi: 3 A j K [ 0 i Aluna l o V 3 I b i*t oh LlNtAR twUATIJNo 

DIMENSION A ( So » :M * UDUh a ( 3 u } » b ( 3 U > 2 ) * U I SO) 

ST AR I 
J - N 

b ( i » 1 ) - A ( J » 2 > 

B ( 1 » 2 ) = A ( L * 3 ) 

D ( 1 ) = C { 1) 

J J = J- 1 
DO b K = 2 » J J 

I h ( ABS ( B ( K. - 1 -> 1 ) ) “Ah S ( 6 ( X ~ i * 2 > ) ) 3 » 4 ♦ 4 
3 B I K > 1 ) = A( K , 2 ) *ti ( N- 1 » 1 ) / I? ( i^-i »<:' ) "ft K t 1 ) 

B l K. 2 ) * A ( R . 3 ) *iM K.~ i » 1 > /o ( N-i i ^ ) 

D ( K ) = ( C ( K ) *b ( • l » 1 ) - a (K . I ) *u ( k - i ) ) / 1> l K, - 1 .2 ) 

GO TO 5 

A B ( K > J ) = A ( . 2 ) - A ( K . 1 ) *0 ( K.-1 » 2 ) / o ( n- t . i) 

B ( K . 2 ) = A { K. . j ) 

D(K)=cns)-Aikin*D(N-n/oiN-i*i^ 
b CONTINUE 

7 X ( J ) - ( C ( J ) * B ( J- I-*n-A(J S il)*UtJ-n)/(A(j.b)‘ot^-l»l)-AlJ.l)*BIU~1.2 
1 ) I 




K -- J ~ 1 

I r ( K — 1 ) 1 0 u » 100, lU 

10 K. = K - 1 

GO TO 15 
TOO RETURN 
END 




^IIT? 
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Card 1 F ormat ( 5F 1 0, 4 , 315 ) 


B1 

- 

Ratio of B^/B. 

ALAM 

- 

Bearing number, 

HI 

- 

Normalized gas film thickness, H^. 

Tol 

- 

Smallness for convergence. 

NL 

- 

Total grid number for the step pad including two 
end points. 

N2 

- 

Total grid number for the left edge of the step pad 
with H « H^. 

Card 2 

Format (8F1G.5) 


VX(X) Array of the squeeze numbers to be solved. 


Card 3 

DX(X) 



Format (8F10.5) 

Array of increments of X defined as 
N , , 

such that Z BX(I) - 1. 

1=1 

Format (8F10.5) 


P0(I) 


Array of the guessed zero order pressure profile Po, i 
on the pad with Po, 1 = Po, n = 1. 



